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I.  INTRODUCTION 

Computer  programs  for  the  new  E-field  solution  for  a  conducting 
surface  of  revolution  are  described  and  listed  in  this  report.  The  new 
E-field  solution  was  introduced  in  [1].  It  is  assumed  that  the  reader 
is  familiar  with  [1].  The  conducting  surface  of  revolution  is  called  S. 
The  surface  S  may  be  either  open  or  closed.  Neither  fins  nor  wires  are 
attached  to  S.  If  S  is  open,  then  S  is  assumed  to  be  bounded  by  at  most 
two  contours,  one  at  the  beginning  of  the  generating  curve  of  S,  and  one 
at  the  end  of  the  generating  curve  of  S. 

[1,  Eq.  (86)]  can  be  rewritten  as 


where 


Z  I  =  V  ,  n  »  0,  ±1,  ±2,... 
n  n  n 


,mm 

zme 

n 

'n 

,em 

zee 

n 

’n 

Equation  (1)  is  called  the  moment  equation,  is  called  the  moment  matrix, 

and  V  is  called  the  excitation  vector.  The  subroutine  ZMAT  of  Section  IT 
n 


M2  where  Ml  and 


calculates  the  elements  of  Z  for  n  =  Ml,  Ml+1,  Ml+2,... 

n 

M2  are  non-negative  integers  and  M2  ^  Ml.  The  subroutine  ZMAT  calls  the 
function  BLOG  of  Section  III.  The  subroutine  PLANE  of  Section  IV  calcu¬ 
lates  the  elements  of  V  for  n  =  Ml,  Ml+1,  Ml+2,... M2.  The  subroutines 

n 

DECOMP  and  SOLVE  of  Section  V  solve  (1)  for  I  .  The  subroutines  ZMAT  and 

n 

PLANE  are  similar  to  the  subroutines  listed  on  pages  51-55  and  61-62  of  [2] 
The  subprograms  BLOG,  DECOMP,  and  SOLVE  are  exactly  the  same  as  in  [2]. 

The  main  program  of  Section  VI  calls  the  subroutines  ZMAT,  PLANE, 
DECOMP,  and  SOLVE  in  order  to  calculate  the  surface  density  of  electric  cur 
rent  and  electric  charge  induced  on  S  by  a  plane  wave  that  propagates  along 
the  z  axis.  The  z  axis  is  the  axis  about  which  the  generating  curve  of  S 
is  rotated.  Such  a  plane  wave  is  called  an  axially  incident  plane  wave. 

For  an  axially  incident  plane  wave,  it  suffices  to  solve  (1)  only  for  n*l. 
Because  the  subroutines  ZMAT  and  PLANE  are  designed  for  n  =  Ml,  Ml+1, 
Ml+2,... M2,  these  subroutines  are  more  general  than  the  main  program  of 
Section  VI. 

In  [1,  Figs.  1,2,  and  3],  the  new  E-field  solution  for  the  current 
and  charge  on  a  small  circular  disk  is  compared  with  Bouwkamp's  power 
series  solution  [3].  In  Section  VII,  Bouwkamp's  power  series  solution  for 
the  surface  density  of  electric  current  on  a  conducting  circular  disk  of 
unit  radius  excited  by  an  axially  incident  plane  wave  is  converted  to  mksc 
units  [A,  p.  1]  for  a  disk  of  arbitrary  radius  a.  The  electric  charge  on 
the  disk  is  extracted  from  this  electric  current.  A  computer  program  which 
calculates  this  current  and  charge  is  described  and  listed  in  Section  VII. 

In  [1,  Figs.  A, 5,  and  6],  the  new  E-field  solution  for  the  current 
and  charge  on  a  small  sphere  is  compared  with  the  Mie  series  solution 


[4,  Eq.  (6-103)],  [5]  for  the  surface  density  of  electric  current  and 
electric  charge  on  a  conducting  sphere  excited  by  a  plane  wave,  A. 
computer  program  which  calculates  the  Mie  series  solution  for  the  cur¬ 
rent  and  charge  is  described  and  listed  in  Section  VIII. 


n , i+1  n ,  j  +1  n  'i+l.j+l 


Ou 


kp  (a**)..  +  a  .(ZH)  -  a  (Z**) .  ) 

i  n  'ij  nj  n  lj  n,j+l  n  i,j+l 


(zme) 

(n  'ij 


kp.((Zt4>)..  -  a  .  (Z*)  .  .  +  a  .  (Z**) .  ,) 

J  n  ij  ni  n  ij  n,x+lv  n  i+l,j 


(zee)..  =  k2p.p.(Z#).. 

n  'ij  KiMjv  n  ij 


where  the  Z^'s  on  the  right-hand  sides  of  (7)— (10)  are  given  by  [2, 
Eqs.  (9)— (12) ]  and  is  given  by  [1,  Eq.  (99)]. 


nj  -  „4.  .  ' 

All  the  subscripts  j  in  (11)  coincide  with  each  other,  but  the  j  which 
multiplies  p^  on  the  right-hand  side  of  (11)  is  /^I.  It  is  evident  from 
[2,  Eqs.  (9)-(12) ]  that 


Equations  (7 )— (12)  imply  that 


Therefore,  it  suffices  to  calculate  the  elements  of  Z^  of  (2)  only  for 
non-negative  values  of  n. 
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New  matrices  (Z^,  n=0,  ±1,  ±2,...}  are  defined  by 


0  BZ, 


where 


Z  =  BZ  ,  n  =  ±1,  ±2,, 
n  n 


B  =  c|  kA1)"2 


Later,  it  will  become  evident  that  the  scale  factor  S  prevents  the 
magnitudes  of  the  elements  of  {Z^,  n  =  0,1,2,...}  from  becoming  exces¬ 
sively  small  as  k  approaches  zero.  Knowledge  of  the  elements  of 

{Z  ,  n  *»  0,1,2, ... }  is  equivalent  to  knowledge  of  the  elements  of 
n 

{Z  ,  n  =  0,  ±1,  ±2,...}  of  (2). 
n 

The  subroutine  ZMAT(M1,M2,NP>NPHI,NT,RR,ZH,X>A.,XT>AT,Z)  is  listed 
at  the  end  of  this  section. .  The  subroutine  ZMAT  puts  the  elements  of  Zn 
of  (14)  and  (15)  in  Z ( (n-Ml)*N*N+l)  through  Z ( (n-Ml+l)*N*N)  for  n  -  Ml, 
Ml+l,Ml+2, . . .M2  where  Ml  and  M2  are  non-negative  integers  and  M2  Ml. 

Here»  n  =  2*NP  -  3  (17) 

Storage  of  Z^  is  by  columns.  Z  is  the  only  input  argument.  The  rest  of 
the  arguments  are  input  arguments  and  have  the  same  meanings  as  in  the 
subroutine  ZMAT  listed  in  [2,  pp.  51-55].  NP  is  the  number  of  data  points 
on  the  generating  curve  of  S.  RH  and  ZH  contain  the  coordinates  of  the 
data  points  in  electrical  length. 


J=1,2,...NP 

(18) 

J=1,2,...NP 

19) 

Here,  k  is  the  propagation  constant,  z  is  the  coordinate  along  the  axis 
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about  which  the  generating  curve  of  S  is  rotated,  and  p  is  the  distance 
from  this  axis.  Also,  (t^)  denotes  evaluation  at  the  Jth  data  point. 
NPHI  is  n.  in  the  Gaussian  quadrature  formulas  [2,  Eqs.  (64 ) — (66) ] .  X  , 

(nJ 

contains  the  n^  abscissas  ^  of  [2,  Eq.  (70)],  and  A  contains  the 
(n. ) 

n.  weights  A.  ™  in  [2,  Eqs.  (64)-(66)].  NT  is  n  in  the  Gaussian 

(f>  >C  t 

quadrature  formulas  [2,  Eqs.  (62)-t63)].  XT  contains  the  n^  abscissas 

(nt)  (nt) 

x^,  in  [2,  Eqs.  (62)-(63)],  and  AT  contains  the  nt  weights  A^,  in 

[2,  Eqs.  (62 )— (63 ) ] -  The  subroutine  ZMAT  calls  the  function  BLOG  which 

is  listed  in  the  next  Section. 

Minimum  allocations  in  the  subroutine  ZMAT  are  given  by 


COMPLEX  Z (M3*N*N) ,  G4A(M3),  G5A(M3),  G6A(M3), 
G4B(M3) ,  G5B(M3),  G6B(M3),  GA(NPHI) ,  GB(NPHI) 


DIMENSION  RH(NP) ,  ZH(NP),  X(NPHI),  A(NPHI), 

XT (NT),  AT (NT),  RS(NP-l),  ZS(NP-l),  D(NP-l), 
DR(NP-l) ,  DZ(NP-l),  DM(NP-l) ,  C2(NPHI), 

C3 (NPHI) ,  R2 (NT) ,  Z2(NT),  C4(M3*NPHI),  C5(M3*NPHI), 
C6(M3*NPHI) ,  Z7 (NT) ,  R7 (NT) ,  Z8(NT),  R8 (NT) 


where  N  is  given  by  (17)  and 

M3  =  M2  -  Ml  +  1  (20) 

In  view  of  (5),  (14)  becomes 

i*  °  ' 

(21) 

o 

where  the  ijth  elements  of  Z^  and  Z ^  are  given  by  [2,  Eqs.  (9)  and  (12)]. 


B 


In  [2,  Eqs.  (9)  and  (12)],  the  region  of  integration  for  which 

S+i 


t“  <  t  <  t  ^  I 
p  -  -  P+1  [ 


(22) 


t  <  t’  <  t  , . 
q  -  -  q+1 


is  called  A  .  The  ranges  of  values  of  the  subscripts  p  and  q  on  A  are 


P9 


pq 


given  by 


(23) 


I  <  p  <  p-l 

I I  q  1  p-i  J 

where,  as  in  [2],  P  is  the  number  of  data  points  on  the  generating  curve 
of  S.  From  [2,  Eq.  (48)],  the  contributions  to  the  elements  of  Z^*"  due  to 
A  are 

pq 

*  3^  A 

(Z„), .  =  — — P  q  (G_  sin  v  sin  v  +  G_  cos  v  cos  v  )  + 

0  ij  8  5a  p  q  7a  p  q7 


(-l)q"j  jk2A  A 


2_3 


(Gr.  sin  v  sin  v  +  G_,  cos  v  cos  v  )  - 
u  5b  p  q  7b  p  q 

(_1)p+q-l-j  i  (24) 

where  i  is  either  p-l  or  p  and  j  is  either  q-1  or  q  but  neither  i  nor  j 

can  be  0  or  P-l.  The  asterisk  on  the  left-hand  side  of  (24)  denotes 

that  (z!it),.  is  not  all  of  (Z^1")  .  ,  but  only  the  contribution  due  to  A 

0  ij  0  ij  pq 

According  to  [2,  Eq.  (51)],  the  contribution  to  the  elements  of  32^  due 


to  A  is 

pq 


2 

,  ,  8k  A  A  A  sin  v 

<\q  ■  «  <-l^,(C5a  +  G5b> 


44  > 


(25) 


Equation  (25)  means  that  A  contributes  only  to  3(ZI/ )  and  gives  all 

pq  U  PQ 


°£ 


pq 


y  i 


8 


Substitution  of  (2)  into  (15)  gives 


Z  = 
n 


BZ™ 

n 

Bzr 

Sz 

n 

8znee 

„mm 

„em  „me 

z  , 

Z  ,  Z 

n 

n  ’  n 

,  n=l,2,3, . . . 


(26) 


;ee 

n 


(9),  and  10,  respectively.  The  contributions  to  the  elements  of  BZ1™1  due 

n 

to  A  are 
pq 

B(zmm)  =  6((ztt)  -o  (Zt4>)  +  ol  -a  a  (Z^)  ) 

n  p-l,q-l  n  p-l,q-l  nq  n  p-l,q  np  n  p,q-l  np  nq  n  pq 

(27) 

8(Zmni)  .  =  B((ZCt)  .-a  (Zt4>)  -  a  (Z^)  +  a  a  (Z^)  )  (28) 

n  p,q-l  n  p,q-l  nq  n  pq  np  n  p,q-l  np  nq  n  pq 

B(^)  1  fl  =  6<(*nt)n  1  a  +  1  n  +  +  <*  <*  (Z!<t>)na)  (29) 

n  p-l,q  n  p-l,q  nq  n  p-l,q  np  n  pq  np  nq  n  pq 


ea""")  -  B((Ztt:)  +  a  (Zt4>)  -a  (Z*C)  -a  a  (ZH)  ) 

n  pq  n  pq  nq  n  pq  np  n  pq  np  nq  n  pq' 


(30) 


The  asterisks  above  the  Z  's  in  (27)- (30)  denote  the  contributions  due  to 

n 

eiii 

A  .  The  contributions  to  the  elements  of  8Z  due  to  A  are 
pq  n  pq 


8 (*em)  .  =  8kp  ((*<!>t)  .  «  , 

n  p,q-l  p  n  p,q-l  nq  n  pq 


-  a  _) 


(31) 


8(Zem)  =  8kp  ((Z^)  +«nn(z^>nn) 

n  pq  p  n  pq  nq  n  pq 


(32) 


me 

The  contributions  to  the  elements  of  BZ  due  to  A _  are 

n  pq 


8(Zme)  .  -  8kp  ((Zt<J>)  .  +  a  (Z^)  ) 

n  p-l,q  q  n  p-l,q  np  n  pq 


(33) 


8(Zme)  =  8kp  ((Zt<}>)  -  a  (Z^)  ) 

n  pq  q  n  'pq  np  n  'pq' 


(34) 
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The  asterisks  above  the  Z  's  in  (3 1) — (34)  denote  the  contributions  due 

n 

66 

to  A  .  As  concerns  Z  ,  A  contributes  only  to  the  path  element  of 
pq  n  pq 

66 

Z  and  gives  all  of  this  element. 

n 


8(Zee)  =  Sk2p  p  (z*> 
n  'pq  Mp^q  n  pq 


In  (23),  the  subscripts  p  and  q  on  A  run  from  1  to  P-1.  However, 

i=l,2,...P-2  on  the  testing  function  of  [1,  Eq.  (102)],  and  j=l,2,...P-2 

on  the  expansion  function  Jm.  of  [1,  Eq.  (100)].  Therefore,  some  of  (27)- 
(34)  must  be  deleted  when  p  is  either  1  or  P-1  or  when  q  is  either  1  or  P-1. 

If  p=l,  then  (27), (29),  and  (33)  are  to  be  deleted.  If  p  =  P-1,  then  (28), 

(30), and  (34)  are  to  be  deleted.  If  q*l,  then  (27),  (28),  and  (31)  are  to 
be  deleted.  If  q  =  P-1,  then  (29),  (30),  and  (32)  are  to  be  deleted. 

The  Z^'s  on  the  right-hand  sides  of  (27)— (35)  are  given  by  [2,  Eqs. 
(19 )— (22) ] .  If  [2,  Eqs.  (19)— (22)]  are  substituted  into  (27)-(34),  then, 
thanks  to  (11)  and  the  formulas 


_d_ 

dt 


Tp(t) 


P  (t) 


_d_ 

dt 


(36) 

(37) 


the  last  term  in  [2,  Eq.  (19)]  and  the  terms  in  [2,  Eqs.  (20)— (22) 1 
cancel  each  other.  These  terms  are  the  electric  scalar  potential  terms. 

The  remaining  terms  in  [2,  Eqs.  (19) - (22) ]  are  the  magnetic  vector  potential 
terms.  Therefore,  [2,  Eqs.  (19) — (22) ]  reduce  (27) — (34)  to 


B(Z  )  .  .  =  B((Z  )  ,  ,  -  a  (Zwr)  .  +  a  (ZYl-)  ,  -  a  a  (ZTr)  ) 

n  p-l,q-l  n  p-l,q-l  nq  n  p-l,q  np  n  p,q-l  np  nq  n  pq 


BCZ™)  =  6((ZCt)  -a  (Zt4>)  -a  (Z^)  +  a  a  (Z^)  )  (39) 

n  p,q-l  n  p,q-l  nq  n  pq  np  n  p,q-l  np  nq  n  pq 

BCZ™1)  .  =  3  ( (ZtC)  +  a  (Ztt(>)  +  a  (z^*)  +  a  a  (Z^)  )  (40) 

n  p-1  ,q  n  p-l,q  nq  n  p-l,q  np  n  pq  np  nq  n  pq 


B(Zmm)  =  6((Ztt:)  +  a  (Z*^)  -a  d^)  -a  a  (Z^)  ) 

n  pq  n  pq  nq  n  pq  np  n  pq  np  nq  n  pq 


B(z“)  =  Bkp  n((zf)  .-a  (!J*)  ) 

n  p,q-l  p  n  p,q-l  nq  n  pq 


8(2®")  =  Bkp  D((Z^t)  +  a  (zf)  ) 

n  pq  p  n  pq  nq  n  pq 


8(Zme)  .  =  Bkp  ((!*♦)  +  a  (!*»)  ) 

n  p-1 ,q  q  n  p-l,q  np  n  pq 


B(Z™e)  =  ekp((Z^)  -  a  (Z^)  >  (4! 

n  pq  q  n  pq  np  n  pq 

a 

where  the  Z  's  are  the  magnetic  vector  potential  contributions  to  the  Z  's 
n  n 

of  [2,  Eqs .  (19)-(22)].  Equation  (35)  remains  unchanged. 

B(Zee)  =  Bk2p  p  (ZW)  (46 

n  'pq  "pq  n  pq 

As  extracted  from  [2,  Eqs.  (19) — ( 22) ] ,  the  Z^'s  on  the  right-hand 
sides  of  (38)-(46)  are 


<!“>. 


rp+1  rq+i  2 
i  dt  dt'  kA 


dt'  k  TJ (t)T^ (t1) (Gcsin  v  sin  v' +G,cos  v  cos  v') 
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(!*>  = .  ± 

(  n  }pj  Pp  J 


P+1 


dt  P  (t) 


P  J- 


q+1 


dt'(k2p  T.(t')  G,sin  v') 


(48) 


(ft*)  =  J_ 

(  n  iq  Pq  J 


fP+1 


dt 


q+1 


dt'  P  (t ' ) (k2p'  T . (t)  G,sin  v) 
q  10 


t  t 

p  q 


(49) 


(g<Ms  _  j 

n  pq  PpPq  J 


P+1 


dt  Pp(t)  j  dt'  Pq(f)(k2pp'  g5) 


q+1 


(50) 


(Z^)  =  -J- 


fcp+i  rq+i 


n  Pq  PpPq  j 


dt  P  (t)  dt'  P  (t')(k2pp'  G,  -  n2G_) 

p  q  57 

t  t 

p  q 


(51) 


where  i  is  either  p-1  or  p  and  j  is  either  q-1  or  q,  but  neither  i  nor  j 
can  be  0  or  P-1. 

If  the  approximations  that  led  from  12,  Eqs.  (19)— (22) ]  to  [2,  Eqs. 
(48)— (51) ]  are  applied  to  (47)-(51),  then  (47)-(50)  reduce  to  the  vector 
potential  terms  in  [2,  Eqs.  (48)— (51) ]  and  (51)  reduces  to  [2,  Eq.  (51)]. 
Hence , 
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n  pq  n  'pq  J  2p  2p  7a 


(55) 


(56) 


where  i  is  either  p-1  or  p  and  j  is  either  q-1  or  q,  but  neither  i  nor  j 
can  be  0  or  P-1. 

The  developments  that  have  occurred  from  (21)  to  (56)  are  summarized 

as  follows.  It  was  shown  that  is  given  by  (21)  where  the  contributions 

to  the  elements  of  Z^C  due  to  A  are  given  by  (24)  and  the  elements  of  gZ^ 

0  pq  0 

are  given  by  (25).  Equation  (26)  was  established  for  2^.  The  contributions 
to  the  matrix  elements  on  the  right-hand  side  of  (26)  due  to  are  given 
by  (38)-(46).  The  zn's  on  the  right-hand  sides  of  (38)-(46)  are  given 
by  (52)- (56). 

The  subroutine  ZMAT  is  similar  to  the  subroutine  described  and  listed 

in  [2,  pp.  43-55].  Hence,  it  suffices  to  point  out  statements  that  differ 

from  those  in  the  subroutine  listed  in  (2,  pp.  51-55].  In  the  listing  of 

ZMAT  at  the  end  of  this  Section,  line  53  puts  3  of  (16)  in  RD. 

The  index  JQ  of  DO  loop  15  obtains  the  subscript  q  on  A^  of  (22). 

This  q  appears  on  the  right-hand  sides  of  (24)— (25) ,  (38)-(46),  and  (52)— (56) . 

With  reference  to  (44)  and  (45) ,  line  59  puts  3kp^  in  RQ.  The  index  IP  of 

DO  loop  16  obtains  the  subscript  p  which  appears  on  the  right-hand  sides 

of  (24)— (25) ,  (38)-(46) ,  and  (52)-(56).  With  reference  to  (42)  and  (43), 

2 

line  81  puts  gkp^  in  RP.  Line  260  puts  gk  p^p^  of  (46)  in  RPQ. 

Inside  DO  loop  31,  the  elements  of  7,  are  obtained  for 


n  =  Ml-1  +  M 


(57) 


where  M  is  the  index  of  DO  loop  31.  Table  1  describes  the  action  of 
some  statements  in  DO  loop  31.  The  statement  whose  line  number  is  given 
in  the  third  column  of  Table  1  stores  the  text  quantity  of  the  second 
column  in  the  variable  in  ZMAT  listed  in  the  first  column.  Integers 
{KM,  M=l,2,...8}  are  defined  by  lines  272-279.  Lines  317-343  set  Z(KM) 
equal  to  VM  for  M=4,6,  and  8,  and  add  VM  to  Z(KM)  for  M=l,2,3,5,  and  7. 
The  branch  statements  among  lines  317-343  prevent  alteration  of  Z(KM) 
at  the  forbidden  values  of  p  and  q  in  [2,  Table  2  on  p.  50]. 
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LISTING  OF  THE  SUBROUTINE  ZMAT 

THE  SUBROUTINE  Z MAT  CALLS  THE  FUNCTION  BLOG 

SUBROUTINE  ZMATI N1 . M2.NP.NPHI.NT.RH. ZH. X.A.XT.AT.Z) 

COMPLEX  ZIl 600) . LI. U2.U3.U4.U5. U6.U7.U8.UA.UB.G4AC 10) .G5A<10) 
COMPLEX  CMPLX*G6A( 10) »G4B< 10).G5B< 10) .G6BC 1 0 ) . H4 A .H5A . H6A ,H4B .MSB 
COMPLEX  H6B.UC.U0.GA( 48) . GBC48 ) • AP. AO, V 1 • V2 • V3 . V4. V5. V6. V7. V8, V9 
01 MENS l ON  RH ( 43 ) • ZH<  43) «X(48).A( 48 ) • XT ( 10).  AT  C 10).RS( 4  2).ZS( 42) 

O I  MENS  I ON  D<42).0R(42).DZ<42).DM(42).C2(48) .C3<48).R2C 10) ,Z2( 10) 
DIMENSION  C41200) .C5{ 200) ,C6( 200 ).Z7( 1 0) .R7« 10) .Z8I 10) .R8( 10) 

C  T =2  • 

CP=.l 

00  10  1=2. NP 
12=1-1 

RSI  I2)=.5*<RH<  I  )  *RH{  12)) 

ZSI  12  )=.5*CZH<  I )*ZHC  12)  ) 

D1  =  .5*(RH{ I  )-RH(  12) ) 

02= .5*1 ZH( I )-ZH< 12) ) 

0(  1 2 )=SQRT( 01*01*02*02) 

OR  I  12 )  =  01 
OZ( 12 )=02 

OM( 12 )=0( I 2 )/RS( 12) 

10  CONTINUE 
M3=M2-M1*1 
M4=  Ml— 1 

P 12=1 .570 796 
OC  11  K=l *  NPHI 
PH=PI 2* (X<K)*1. ) 

C2(K)=PH*PH 
SN=SIN( .5*PH) 

C3(K)=4.*SN*SN 
A  1=P I  2*  A<K ) 

04= .5* A  1 *C3 { K ) 

05= A 1 *COS( PH ) 

D6=A1 *S IN I PH) 

M5=K 

OC  29  M=1.M3 
PHM=(M4*M)*PH 
A 2= COST  PHM ) 

C4<  M5 )  =  D4*A2 
C5IM5  )=05*A2 

C6(M5)=Ob*SlN«PHM)  , 

M5=M5*NPHI 
29  CONTINUE 

11  CONTINUE 
MP=NP-1 
MT=MP-1 
N=M  T*MP 
N2N=MT*N 
N2=N*N 
U1=(0.. .5) 

U2= (0 • . 2. ) 

JN=-1-N 

RD=1. /( D( 1 ) *D(  1 ) > 

00  IS  JQ* 1 »MP 
K0=2 

IF(JQ.EQ.l)  KQ= 1 
1FCJO.EQ.MP)  KQ=3 
R1«RS( JO) 

RQsRl *R0 
Zl-ZS(JO) 


01=01  JO  ) 

D2=0R<  JO) 

03-OZ ( J  Cl 
D4=D2/RI 
05=0M< JQ) 

SV=02/D1 
C V=D3/01 
T6=CT*Ot 
T62=T6*01 
T62=T62*T62 
RC- CP  *R1 
R62=R6*R6 
OG  12  L=1 *NT 
R2IL>=Rl*D2*XTfl_ I 
Z2(L)=Zl+03*XT(t_l 
12  CCNTINOE 
U3=D2*U1 
04=03*01 
00  16  1  P=  1  •  PP 
R3=RS( IP1 
RP=RJ*RO 
23=ZS< IP I 
R4=R1-R3 
24=21-23 
FM=R4*SV*Z4*CV 
PHM=ABS(FM) 

PH=ABS< R4*CV-24*SV1 
OC=PH 

IFfPHM.LE.Dll  GO  TO  26 
06=FHM— 01 

06=SQPTCD6*06+PH*PH) 

26  IF(lP.Ea.JO)  GC  TC  27 
KP=1 

1F(T6.GT.06>  KP=2 
IFIR6.GT.D6)  KP=3 
GO  TO  28 

27  KPx4 

28  GG  TO  141 .42.41, 421.KP 
42  00  40  1=1 «NT 

07=R2(L »-R3 
08=22<Ll-23 
Z7<L 1=0 7*07*08*08 
R7fL)=R3*R2«L) 
Z8(L)=.25*Z7(L) 
R8IL)=.25*R7«LI 
40  CONTINUE 

Z4=R4*R4*24*24 
R4=R3*R 1 
R  5= • 5*R3*SV 
00  33  K=1*NPHI 
A1=C3(K> 

RR=Z4*R4*A1 

UA*0. 

ue=o. 

IFCRR .LT.T62)  GO  TO  34 
00  35  L  =  1  •  NT 
R*SCRT( 27<L)*R7CLl*Al I 
SN=— S  IMP) 

CS»COS(R) 

0C= AT (L )/R*CPPLX (CS»  SN ) 


121 

UA=UA4UC 

122 

UBsXTIL )*OC*UB 

123 

35 

CONTINUE 

124 

GO  TO  36 

125 

34 

00  37  L= 1 .NT 

126 

RsSCRT < Z8(L) 4R8IL J4A1 ) 

127 

SN=-SIN<R| 

120 

CS=COS( R) 

129 

UC=AT<L)/R*SN4CMR_XI-SN,CS 

130 

UA=UA4UC 

131 

UB=XT(L J4UC4UB 

132 

37 

CONTINUE 

133 

A2=FM4R54A1 

134 

D9=RR— A24 A2 

135 

R=AES I A2) 

136 

07=R-D1 

137 

08=9401 

130 

06=  SORT < 08408409 ) 

139 

RsSCRT ( 07407409 ) 

140 

IFID7.GE.0.)  GO  TO  38 

141 

A  1  =  ALCGI ( 08406 ) 4 1— 074R )/D9 

142 

GO  TO  39 

143 

38 

A 1= ALOG ( ( 08406)/ ( 074R ) ) /D1 

144 

39 

UA&A14UA 

145 

U8=A24C 4a/( 064RI— 41 1/014UB 

146 

36 

GA<K)=UA 

147 

GB( K»=UB 

140 

33 

CONTINUE 

149 

•K  1=0 

ISO 

OO  45  M= 1 • M3 

151 

H4A=0« 

152 

H5A=0. 

153 

H6A=0 • 

154 

M4B=0 • 

155 

H5B=0 • 

156 

H6B=0. 

157 

OO  46  K=1.NPHI 

150 

K1=KI41 

159 

D6=C4 ( K 1 9 

160 

07=051 K 1 1 

161 

08=C6( K1 ) 

162 

UA=GA«K» 

163 

UB=GB«K» 

164 

H4A=064UA4H4A 

165 

HSA=074UA4H5A 

166 

H6A=084UA4H6A 

167 

H4B=064UB4H4B 

160 

H5B=074UB4H58 

169 

H6B=084UB4H6B 

170 

46 

CONTINUE 

171 

G4A(MI=H4A 

172 

G5A(M)=HSA 

173 

G6A(MJ=H6A 

174 

G4BIM |=H4B 

175 

G5E(M)=HSB 

176 

G6BC M)*H6B 

177 

45 

CONTINUE 

170 

IFIKP.NE.4)  GO  TO  47 

179 

A2=0!/(PC24R1I 

100 

06=2. /Ol 

20 


u. 

J 


Lt 

lai 

08=0. 

182 

DO  63  Ksl.NPHI 

183 

A1=R4*C2(K> 

184 

R=R4*C3|K) 

1 

185 

IF IR.LT .162 )  CO  TO  64 

r 

186 

07*0. 

t* 

L'- 

187 

00  65  L=1.NT 

' 

188 

07=074  AT  <L I /SORT  f  Z7(L)4A1 | 

189 

65 

CONTINUE 

190 

GO  TO  66 

I 

191 

64 

A1=A2/(X<K)41.) 

V 

192 

07=06* ALOG 1  A 1 4SQRT <1.4At*Al»> 

193 

66 

08=08*A(K)*07 

194 

63 

CONTINUE 

195 

A  1=  «5*A2 

196 

A 2= 1 ./A  1 

i 

197 

08=-P12*0842./Rl*(BL0G« A2)4A2*8L0GCAt>> 

198 

00  67  M= 1 . M3 

199 

G5AC  M  >=084G5Af Ml 

:- 

200 

67 

CONTINUE 

201 

GC  TO  47 

202 

41 

00  25  M=1.N3 

i 

203 

G4A(M)=0. 

j 

204 

G5A(H)=0. 

205 

G6A (M )=0. 

206 

G4B(M)=0. 

207 

G58(M )=0 • 

« 

208 

G6E(M)=0. 

i 

209 

25 

CONTINUE 

I 

210 

00  13  L=1.NT 

211 

A1=R2CL> 

212 

R4=A1— R3 

-* 

213 

Z4=Z2(L)-Z3 

214 

24=P4*R44Z4*Z4 

i 

215 

R4=R3*A 1 

216 

00  17  K=1.NPHI 

217 

RsSCRTl Z44R4*C3<K) ) 

■ 

218 

SN=— SINIRI 

* 

219 

CS=C0S(R» 

220 

GA<K)=CMPLX(CS.SN)/R 

«i 

2?1 

17 

CONTINUE 

222 

06=0. 

223 

IF<  P62.LE •  Z4J  GO  TO  51 

224 

00  62  K=1.NPHI 

225 

06=064 A ( K l/SCRT ( Z44R4*C2(K>1 

Mj 

226 

62 

CONTINUE 

227 

Z  4=  3 • 141 593 /SORT { Z4/R4I 

228 

D6=-P 1 2  *D64 ACOG( Z44SQRTC 1 • 4Z4*Z4) ) /SORT CR4 ) 

229 

51 

A1«AT<L) 

230 

A2=XT«->«A1 

231 

K  1*0 

232 

00  30  M«1.N3 

M  233 

US*0. 

; 

234 

U6=0. 

, 

*35 

U7*0. 

*- 

236 

00  32  K=1.NPHI 

237 

UAsGACK) 

h  938 

K  1*K1 4 1 

M  ?39 

U5=C4(Kl I4UA4U5 

240 

U6=CSCK1»*UA4U6 

n 


r 


ri 

■j 


> 


I 


L 


241 

U7=C6(K1 14UA4U7 

242 

32 

CONTINUE 

243 

U6=D6+U6 

244 

64 A ( M  )=  A1 4U5+G4A (Ml 

245 

G5A(M|=A1 *U6*G5A  <  M ) 

246 

G6A (M)sAl 4U74G6A (Ml 

247 

GAB ( N |= A24U5+G4B 1 Ml 

248 

G5B( M *= A2*U6 ♦G5B < M ) 

249 

G6B(M)=A2*U74G68<M> 

250 

30 

CONTINUE 

251 

13 

CONTINUE 

252 

47 

A1=CRIIP) 

253 

UA=Ai  *U3 

254 

UB=0Z IIP) *U4 

255 

A  2=0( IP) 

256 

06= -A  2*02 

257 

07=01 *A1 

2S8 

08=01 *A2 

259 

JM=  JN 

260 

RPQ=R1*RP 

261 

00  31  N=1 #M3 

262 

FM=M4+M 

263 

HSA=G5A(M» 

264 

H5B=G5Bf M) 

265 

H A A=GAA ( M) 4H5A 

266 

H4B=G4B|M) AH58 

267 

U7=UA*H5A+UB*HAA 

268 

U8=UA*H58+UB4M4B 

269 

US=U7-U8 

270 

U6=U7+U8 

271 

V9=U2*08*(H5A+04*H58) 

272 

K1=IP*JM 

273 

K2=K1 41 

274 

K3=K1 4N 

275 

K4=K2*N 

276 

K5=K24MT 

277 

K6=K44MT 

278 

K7=K34N2N 

279 

K8=K44N2N 

280 

IF ( FM  »NG  •  0  • )  GO  TO  14 

281 

U7=-U1*H4A 

282 

Vl=U54U7 

283 

V2=U5-U7 

284 

V3=U6— U7 

285 

V4=U6+U7 

286 

V5=0  • 

287 

V6=0. 

288 

V7=0. 

289 

V8=0  • 

290 

ZIK8+MT )=R0*V9 

291 

GO  TO  43 

292 

14 

H6A=G6A(M) 

293 

H6B=G68(M) 

294 

U7=D6«( H6A-H6B) 

295 

Ue=D6*( H6A4H6B ) 

296 

UC*07*( H6A+D4AH6B) 

297 

A1=FM40NCIP) 

298 

FM0=FM*O5 

299 

AP= 1 • /A 1*U1 

300 

AC*C1./FM0)4U1 

22 


UD=A0*V9 
V5=U7-U0 
V6=U8 4U0 
V7=AP*VS 
V8-APAV6 
UD=AQ*UC 
V1=RD*C U5-U04V7I 
V2=RO*( U5-UO-V7I 
V3=RD*4  U6  +  UD+V8) 

V  4=RD*  4 U6+U0- V8 > 

VS=RP*VS 

V6=RP*V6 

U0*AP*V9 

V7=RC*(UC*U0) 

Ve=RQ*  (  UC-UD) 

Z<K8*MT)=RPQ* ( V9-A1*FMD*H4A*U2I 
>  GC  TO  <18,20.19),K0 
)  Z<K6)-V6 

IF( IP.EO.l)  GO  TO  21 
2  {K3)-=Z(K3I  4V3 
Z4K7)=Z(K7)4V7 
IFUP.EQ.MPI  GO  TO  22 
Z ( K4 I =V4 
Z(K8)=V8 
GO  TO  22 
»  Z4K5)=Z4K5)4V5 

IFIIP.EO.l)  GO  TC  23 
Z1K1 J=ZCK| )»V1 
Z4K7|=Z4K7)*V7 
IFUP.EO.MP)  GO  TO  22 
I  ZCK2)-=ZCK2)*V2 
24K£)=V8 
GO  TO  22 
i  Z4K5)=Z<K5)+V5 
Z4K6J=V6 

IF( IP  »EQ. 1 )  GO  TC  24 

Z(Kl)=Z4Kll*VI 

Z<K3)=Z4K3)*V3 

ZIK7)=ZCK7>  *-V7 

IF4IP.E0.MP)  GO  TO  22 

Z4K2)=Z(K2>4V2 

Z€K4)=V4 

Z4K8)=V8 

JV=JM+N2 

CONTINUE 

CONTINUE 

JN= JN  +  N 

CONTINUE 

RETURN 

ENO 
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III.  THE  FUNCTION  BLOG 

The  function  BLOG(x)  calculates  log(x  +]/l  +  x^)  for  x  >_  0. 
The  method  of  calculation  is  described  in  [2,  p.  56]. 


OOIC  LISTING  OF  THE  FUNCTION  8L0G 

002  FUNCTION  8L0GIXJ 

003  IFCX.GT..1)  GO  TO  I 

004  X2=X*X 

005  BLOG=f ( .0754X2-. 1 666667 ) 4X24 1 .) 4X 

006  RETURN 

007  1  BLOG=ALOG{ X+SQRTt 1.4X4XJ ) 

008  RETURN 

009  E  NO 


IV.  THE  SUBROUTINE  PLANE 


The  subroutine  PLANE  calculates  the  elements  of  V  of  (4).  Accord- 

n 

ing  to  [1,  Eqs.  (61)  and  (62)],  the  ith  elements  of  the  column  vectors  V™ 

n 

and  on  the  right-hand  side  of  (4)  are  given  by 


Vm.  =  -  w"1.  •  ElnCds 

m  n  j  — ni  - 

S 


V6 .  =  — —  We.  .  Elncds 

m  n  J  — ni  — 


Expressions  (58)  and  (59)  are  calculated  for  E1  equal  to  the  0-polari.;ed 

0 

plane  wave  electric  field  E°  given  by 
6  t  “jk  4  r 

E  =  Ug  kn  e  (60) 

inc 

and  also  for  E  equal  to  the  ^-polarized  plane  wave  electric  field 


Et  given  by 


*♦-4 

-  ^ 


In  (60)  and  (61), 


k  =  -  k(u  sin  0  +  u  cos  0  ) 

— t  — x  t  — z  t 


^0  =  u^cos  Oj.  (63) 

u*T  =  u  (64) 

~i>  -y 

where  0  is  the  angle  of  incidence  and  where  u  ,  u  and  u  are  unit  vectors 
t  — x  — y  — z 

in  the  x,y,  and  z  directions,  respectively.  Also,  £  is  the  radius  vector 


from  the  origin.  The  origin  must  lie  on  the  axis  about  which  the  generating 


curve  of  S  is  rotated  because  this  axis  is  the  z  axis.  The  electric  fields 


(60)  and  (61)  are  the  same  as  [2,  Eqs.  (108)  and  (109)1. 

If  V™  due  to  E®  is  called  V10®  and  if  V®  due  to  E®  is  called 

n  —  n  n  — 

r\  then,  according  to  (58)  and  (59),  the  ith  elements  of  V01®  and  V®® 
n  n  n 

given  by 


^mQ  _  .L 
ni  n 


W®.  •  E°di 
-ni  - 


ni  n 


We.  •  Eeds 
— ni  — 


If  r  due  to  E^  is  called  and  if  Ve  due  to  E^  is  called  V®^,  then 
n  —  n  n  —  n 

the  ith  elements  of  Vm<^  and  Ve<^  are  given  by 

n  n  ' 


rm4>  = 

1 

ni 

n  Jj 

s 

re4> 

kp. 

l 

ni 

n 

w"1  •  E^ds 

-ni  — 


W®.  •  E^ds 
— ni  — 


For  n=0,  the  testing  functions  and  are  given  by 

[1,  Eqs.  (93)  and  (94)].  Hence, 


^Oi  ^Oi 


kpi^i  "  £  i 


where  and  are  defined  by  [2,  Eqs.  (4)  and  (5)].  Thanks  to  (69) 
and  (70),  each  defined  by  (65)-(68)  is  equal  to  one  of  the  V^'s  de 


fined  by  [2,  Eqs.  (113),  (114),  (116),  and  (117)1.  More  precisely, 


are 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 
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„0  t  . 

H  =  -  u,  k  e 
- 


-jk  •  r 


The  incident  magnetic  field  associated  with  the  <j>-polarized  electric 
field  of  (61)  is  called  and  is  given  by 


Substitution  of  [1,  Eq.  (105)]  for  and  for  H*nC  in  [1,  Eq.  (68)] 


gives 


Vmfl  =  -  *  T . (t )  e“jn^(H6  •  n)ds 
ni  n  I  I  l  —  — 


Substitution  of  [1,  Eq.  (105)]  for  u..  and  for  HinC  in  [1,  Eq.  (68)] 


gives 


V*  -  -Mf  T.  (t)  e'^V  •  n)ds 
ni  n  i  _  _ 


If  ds  is  replaced  by  pd<})  dt,  (78)  and  (79)  become 


Q  L,  fi+2 

vm0  =  _  k  [  dt 

ni  n  I 


I 

i<i>l 


e .  ■rt.-j'** 


d<J>(H  •  n)e 


Vmt  =  -  -  f  dtpT.  (t)  d<j>(H^  •  n)e-^n<^ 

ni  n  _  l  —  ~ 


The  radius  vector  £  that  appears  in  expressions  (76)  and  (77)  for 
H®  and  is  given  by 


r  =  u  p  cos  d)  +  u  p  sin  <b  +  u  z 
—  -x  — y  — z 

From  (62)  and  (82),  we  obtain 


-jk  •  r  =  jk(p  sin  0  cos  (p  +  z  cos  0  ) 
— t  —  t  t 
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The  unit  normal  vector  n  that  appears  in  (80)  and  (81)  is  given  by 


n  =  u^cos  v  cos  +  u^cos  v  sin  <p  -  u^sin  v  (84) 

Equations  (76),  (64),  (84),  and  (83)  lead  to 

a  jk(p  sin  0  cos  <p  +  z  cos  0  ) 

H  •  n  =  -  k  cos  v  sin  <j>  e  (85) 

Similarly,  (77),  (63),  (84),  and  (83)  lead  to 

,  jk(p  sin  9  cos  <p  +  z  cos  0  ) 

H  *n  =  k(cos  O^cos  v  cos  <j>  4-  sin  G^sin  v)e  C  (86) 

0  m0 

Now,  H  *  n  of  (85)  is  odd  in  <p  so  that  V  of  (80)  is  even  in  n.  Moreover, 

—  ni 

•  n  of  (86)  is  even  in  <b  so  that  Vm<!,  of  (81)  is  odd  in  n.  Because 
~  m  ni 

of  [2,  Eq.  (114)]  is  odd  in  n,  of  (74)  is  odd  in  n.  Because  V^j*  of  [2, 

ed) 

Eq.  (117)]  is  even  in  n,  of  (75)  is  even  in  n.  Hence, 


^m0  +m<j> 

-n  -n 


ve0  ve<^ 

-n  -n 


^m0 

n  n 


_vefi  ve<t> 

n  n 


,  n-1 , 2 , . . . 


Therefore,  it  suffices  to  calculate  the  elements  of  the  V^'s  in  (87)  only 
for  positive  values  of  n. 


Matrices  V  are  defined  by 
n 


V  -  0 
n 


v 

0 

0 

0 

CD 

f>e 

si* 

-©- 

-+e0 

v 

n 

ve* 

n 

,  n  -  ±1,  ±2,. 
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where  0  is  given  by  (16).  In  view  of  (6)  and  (73),  the  matrices  of 

(14)— (15)  and  Vn  of  (88)-(89)  allow  (1)  to  be  rewritten  as 


where 


Z  I  =  V  ,  n  =  0,  ±1,  ±2,, 
n  n  n 


0 


0 

•+m(|> 

0 

-|m0 

-♦mi}) 

n 

n 

•Je0 

fe<P 

n 

n 

,  n  =  ±1,  ±2, . 


In  (91)  and  (92),  I™9  is  I™  due  to  (E6,H9),  I®9  is  le  due  to  (E9,H9),  I™1’ 

n  n  —  —  u  n  —  —  n 

is  i"1  due  to  (E^,H^) ,  and  Ie<5>  is  Ie  due  to  (E^.H*). 
n  _  _  n  n  —  — 

The  subroutine  PLANE  (Ml, M2 ,NF,NP,NT,RH,ZH,XT,AT,THR,R)  is  listed 

/v 

at  the  end  of  this  section.  The  subroutine  PLANE  puts  the  elements  of  V_^ 
of  (88)  and  (89)  in  R((K-1)*M3*2*N  +  (n-Ml)*2*N+l)  through  R( (K-1)*M3*2*N  + 
(n-Ml+l)*2*N)  for 

0  =  THR(K)  (93) 

Here 


N  =  2*NP-3  (S 

and  M3  =  M2  -  Ml+1  (S 

Also,  n  =  Ml,  Ml+1,  M1+2...M2  where  Ml  and  M2  are  non-negative  integers 
and  M2  >_  Ml.  Furthermore,  K  =  1,2,...NF.  Storage  of  V_n  in  R  is  by 


columns.  R  is  the  only  output  argument.  The  rest  of  the  arguments  are 
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input  arguments  and  have  the  same  meanings  as  in  the  subroutine  PLANE 
listed  in  [2,  pp.  61-62].  NP  is  the  number  of  data  points  on  the  gener¬ 
ating  curve  of  S.  RH  and  ZH  are  defined  by  (18)  and  (19).  NT  is  n^  in 

the  Gaussian  quadrature  formulas  (2,  Eqs.  (132)  and  (133)].  XT  contains  the 
(n  ) 

n^.  abscissas  in  [2,  Eqs.  (132)- (135) ] ,  and  AT  contains  the  n^,  weights 

(O 

A^  in  [2,  Eqs.  (132)  and  (133)].  THR  is  in  radians. 

Minimum  allocations  in  the  subroutine  PLANE  are  given  by 


COMPLEX  R(NF*M3*2*N) ,  FA(M2+3),  FB(M2+3),  FC(M2+3) 
DIMENSION  RH(NP) ,  ZH(NP),  XT (NT),  AT (NT),  THR(NF) , 
CS(NF),  SN(NF),  R2 (NT) ,  Z2(NT) 


where  N  and  M3  are  given  by  (94)  and  (95),  respectively. 

Thanks  to  (71),  the  ith  elements  of  the  column  vectors  and 


in  (88)  are  given  by 


Ve®  =  V1"® 
Oi  Oi 


<  -  ev|M  (97) 

where  and  are  given  by  [2,  Eqs.  (113)  and  (117)].  If  the  con- 

e0  -  - 

tribution  to  V_.  of  (96)  due  to  the  integration  from  t  to  t  is  called 
Oi  p  p+1 

V®®,  then,  from  [2,  Eq.  (124)], 


.  „  j ir  kA  sin  v  cos  0^ 
See  _  P _ P  t 


ukA  cos  v  sin  6 


<Fla  '  F-la>  ' 


F0a  + 


.  iTrkA  sin  v  cos  0  irkA  cos  v  sin  9 

(.1)P-i( _ E - E - £.  (p  _  F  \ _ 2 - E - 1  v 

'  L>  '  4  'lb  -lb'  2  Ob' 


where  i  is  either  p-1  or  p,  but  i  is  neither  0  nor  P-1.  Here,  P  is  the 


number  of  data  points  on  the  generating  curve  of  S.  The  integration  from 
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t  to  t  , ,  contributes  to  BV„T  of  (97)  only  for  i  =  p  and  gives  all  of 
p  p+1  Oi 

gV^.  From  [2,  Eq.  (127)],  we  obtain 


jgTTkA 


A  sin  v 


- E.  ((p  _  p  )  +  - E.  (p  _  p  ))  (99) 

2  U  la  -la'  2p  K  lb  -lb"  V  ' 

P 


y  P|l(t 

Consider  the  elements  of  BV  and  BV  on  the  right-hand  side  of 

n  n 

(89).  In  view  of  the  integral  formula 


.-n  j  (k  p  sin  0  cos  <j>  +  n<{)) 

J  (k  p  sin  0.)  =  4—  e  d(]> 

n  t  £n 


(100) 


for  the  Bessel  function  J  deduced  from  [6,  Eq.  (9.1.21)],  substitution 

n 

of  (85)  into  (80)  gives 


„  fi.n  2  (i+2 

6  V*  - - —  dtpT  (t)  cos  v  (J  +  J  )e 

m  n  i  n+i  n-i 


jkz  cos  0 


(101) 


where 


J  =  J  (kp  sin  0 ) 
n  n  t 


(102) 


Similarly,  substitution  of  (86)  into  (81)  gives 


.  Q.n+1  .  2  ri+2 

n -  dtpT  (t)  (cos  v  cos  0 

ni  n  J_  i  t 


(Jn+l  '  J„-l) 


-  2  i  sin  v  sin  0  J  )  e 
J  t  n 


jkzcos  0^ 


(103) 


The  contribution  to  BV  .  due  to  the  integration  from  t  to  t  . 

ni  P  P+1 

is  called  B*m(?.  The  contribution  to  BVmf  due  to  the  integration  from 
ni  ni 

t  to  t  is  called  8V™t.  Now, 
p  p-rl  ni 
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p  =  p  +  (t-t  )  sin  v 
P  P  P 


z  =  z  +  (t-t  )  cos  v 
P  P  P 


.  „  .  2 (t-t  ) 

Ti(t)  =  \  [1  +  — a — E  ]  ,i-p-l,PJ 


t  <  t  <  t 
P  ~  ~  P+1 


(104a) 


(104b) 


(104c) 


(104d) 


where  v  ,  p  ,  z  ,  and  t  are,  respectively,  the  values  of  v,  p,  z,  and  t 
P  P  P  P 

midway  between  t~  and  t  Substitution  of  (104)  into  (101)  and  (103) 

P  P+1 


gives 


6jn7Tk2A  p  COS  V 


A  sin  v 


(F  ^  „  +  F  ,  .  + 


(F  u  +  F_  ,  .  ))  - 


n+l,a  n-l,a  2pp  n+l,b  n-l,b 


(-1)P  i6jnTrk2A  p  cos  v 


A  sin  v 


£-= - £  (F  ,,  +F  ,  ))  (105) 


/] p  4_F  4.  _ t.  (V  4-V 

'  n+l,b  n-l,b  2pp  n+l,c  n-l,c 


. n+1  ,  2,  Q 

*  .  gj  7rk  A  p  cos  v  cos  6 
*m<b  PMP  P  t 


A  sin  v 

(v  -F  +  - 2.  (p  -F  ))  + 

'  n+l,r  n-l,a  2pp  1  n+l,b  n-l,b;i 


gin  Fk2A  p  sin  v  sin  0^  A  sin  v  fc-1 

- - E-E - E - E  (F  +  -P  _E  F  ,  +  -H 

2n  2p  nb  -3 

' 

gin+*TTk^A  p  cos  v  cos  0  A  sin  v  ;.vff 

/  nP-i  /_ _ P  P  P _ t  fp  _F  +  _2 - E  (F  -F  ))  + 

'  '  4n  '  n+l,b  n-l,b  2Pp  n+l,c  n-l,c^  *3 


0jn7Tk2A  p  sin  v  sin  G 
P  P  P  t 


A  sin  v 

(F  +  -B— - E.  F  )) 

nb  2P  nc 

P 


(106) 
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where 


„  fP+1  jkz  cos  0 

-7—  J  (k  p  sin  0  )  e  C  dt 

A  m  t 


jkz  cos  0 


(t — )  (t-t  )  J_(kpsin  0  )  e 

A  J  _  pm  t 

F  t 

P 


(107a) 


(107b) 


„  -  rp+1  „  jkz  cos  0 

(-7—)  (t-t  )  J  (k  p  sin  0  )  e  dt 

A  pm  t 

p 

p 


(107c) 


In  (107),  m  =  n-l,n,n+l.  Also,  p  and  z  are  given  by  (104b)  and  (104c). 

Thanks  to  (87),  expressions  (105)  and  (106)  can  be  viewed  as  8V  .  and 

-ni 

PV-ni’  resPectively* 

Evaluation  of  (107)  by  means  of  an  n^,  point  Gaussian  quadrature 
formula  yields 


"l  (O 


Fma  '  AH  Vk  <Vln  0t>  e 


jkz^  cos  0t 


(108a) 


T  (n  )  (n  )  jkz„  cos  0 

Fmb  =  J=1  kl  x£  Jm(k  P£stn  V  e 


(108b) 


T  (iO  (n  )  ?  jkZp  cos  0 

Fmc  =  A£  (x£  >  Jm(k  P£sin  V  e  (108c) 

where,  as  in  [2,  Eqs.  (134)  and  (135)], 

(n  ) 

A  x5 

p  =  p  +  -&—Z -  sin  v  (109) 

v  p  2  p 

(n  ) 

Vs. 

z„  =  z  +  -  ■z -  cos  v  (110) 

S  p  2  p 

(nT)  (nT) 

The  abscissas  x0  and  weights  A„  in  (108)- (110)  are  the  same  as  in 


[2,  Eqs.  (132) — (135) ] . 


I 


34 


Thanks  to  (74)  and  (75),  the  elements  of  the  column  vectors  $\T 

—►"0(f) 

and  BVn  on  the  right-hand  side  of  (89)  are  given  by 


3ve®  =  Bkp. 

ni  l  ni 


(111) 


sve?  -  Bkp.  V*? 

ni  l  ni 


(112) 


where  and  are  given  by  [2,  Eqs.  (114)  and  (117)].  The  integration 

from  t  to  t  , ,  contributes  to  only  for  i=p  and  gives  all  of  Upon 

p  p+1  ni  np 

replacement  of  i  by  p  in  (111)  and  substitution  of  [2,  Eq.  (125)]  for  V^, 

np 

the  negative  of  (111)  becomes 


BjnTtk2A  p  cos  0 


A  sin  v 


<F  .  i  +F  ,  + 

n+1,  a  n-l,a 


(F„+i,b+F„-i,b»  <113) 


where  F  and  F  ^  are  given  by  (108a)  and  (108b).  Similarly,  substitution 
of  [2,  Eq.  (127)]  into  (112)  yields 


„ . n+1  .  2 . 

Bj  TTk  A  p 


(F  -F  ,  + 

n+1, a  n-l,a 


A  sin  v 

_E — _ — E  (p  _f  )) 

2p  k  n+l,b  n-l,b" 
P 


(114) 


Thanks  to  (87),  expressions  (113)  and  (114)  can  be  viewed  as  BV_np  and 
Bvf^p ,  respectively. 

In  the  subroutine  PLANE,  the  elements  of  of  (88)  are  calculated 
by  means  of  (98)  and  (99),  and  the  elements  of  of  (89)  are  calculated  for 
negative  values  of  n  by  means  of  (105),  (106),  (113),  and  (114).  The  index 
IP  of  DO  loop  12  obtains  p  in  (98),  (99),  (105),  (106),  (113),  and  (114). 

DO  loop  13  puts  ^  pf  of  (109)  and  kz^  of  (110)  in  R2(L)  and  Z2(L),  respec¬ 
tively,  for  P.  =  L.  The  index  K  of  DO  loop  14  obtains  the  Kth  angle  of 


incidence  0  of  (93). 
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The  index  L  of  DO  loop  15  obtains  A  in  (108).  Lines  57-82  calcu¬ 
late  S  and  BJ(m+2)  so  that 


BJ (m+2 )  =  S*J  (ko„  sin  0  ),  m  =  Ml-l.Ml,  ...  M2+1 
m  x-  t 

m  ^  -1 

The  calculation  of  BJ(m+2)  and  S  is  described  in  [2,  p.  59].  As  the  in¬ 


(115) 


dex  L  of  DO  loop  15  changes,  line  88  accumulates  F  of  (108a)  in  FA (m+2), 

m3 

line  89  accumulates  F  ,  of  (108b)  in  FB(m+2),  and  line  90  accumulates  F 

mb  me 

of  (108c)  in  FC(m+2).  If  F_^,  F  ,  and  F_^  are  needed,  lines  94-96  use 
the  formulas 


F  ,  =  -  F. 

-la  la 


F-lb  =  "  Flb 
F-lc  =  "  Flc 


(116) 


to  store  F_^a,  F_^»  and  F_^c  in  FA(1),  FB(1),  and  FC(1)  respectively. 

In  DO  loop  27,  (98)  and  (99)  are  obtained  if  M  is  2.  If  M  is  greater 
than  2,  then  (105),  (106),  (113),  and  (114)  are  obtained  for  n=M-2.  If  M  is 
2,  then  R(K1),  R(K2),  R(K2+MT),  R(K4),  R(K5),  and  R(K5+MT)  are  reserved  for 

tfcgQ  &00  md) 

V0  p_^>  V0p,  0,  0,  0,  and  |3V  ,  respectively.  If  M  is  greater  than  2,  then 

R(K1) ,  R(K2) ,  R (K2+MT) ,  R(K4),  R(K5),  and  R (K5+NT)  are  reserved  for  8Vm6  ,, 

n )  p-1 

8Vm9,  -8Ve0,  -3Vm^  ,  -3Vm<^,  and  ,  respectively.  Table  2  describes  the 

np  np  n,p-l  np  np 

action  of  some  statements  in  DO  loop  27.  The  statement  whose  line  number  is 
given  in  the  third  column  of  Table  2  stores  tne  text  quantity  of  the  second 
column  in  the  variable  in  PLANE  listed  in  the  first  column. 
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Table  2.  Variables  in  PLANE  versus  text  quantities. 
Variable  Text  Quantity 


Variable 
in  PLANE 


F  -  F 

n+l,a  n-l,a 

F  -  F 

n+1 ,b  n-l,b 

.n+1 
J  tt 

j-nrkA  sin  v  cos  0 


Number 


101 

102 


irkA  cos  v  sin  0 


jTrkA  sin  v  cos  0 


<Fla-F-la>  - 


TTkA  cos  v  sin  0 
P  P  « 
2 


jirkA  sin  v  cos  0  TikA  cos  v  sin  0 

- _ E _ E _ £  (F  _F  v _ _2 _ 2 _ t 

4  1  lb  -lb'  2 


KK5+MT)  6  V 

F2A  F 


< 

F  , 

+ 

F  , 

n+1,  a 

n-1  ,a 

F 

+ 

F  ,  . 

n+1  ,b 

n-1  ,b 

F  , 

+ 

F  , 

n+l,c 

n-1  ,c 

F 

F 

n+l,c 

n-l,c 

0jnrk^A  p  cos  v 


B.inTik^A  p  cos  v 


BjnTrk^A  p  cos  v 


A  sin  v  I 

-(F  +F  +  - Ef  f  +F  )) 

1  n+1, a n-1, a  2p  n+l,b  n-l,b  " 

P 


A  sin  v 

(F  +F  +  -2- - E(F  +F  )) 

'  n+l,b  n-l,b  2pp  v  n+l,c  n-l,c" 


R(K2+MT) 

R(K5+MT) 


38 


:  ooic 

LISTING  OF  THE  SUBROUTINE  PLANE 

i-'  002 

SUBROUTINE  PLANE  C  Ml «  M2 • NF  «  KP  • NT • RH 

003 

COMPLEX  RI24C ).U*U1 .UA.UB ,UC .FA f 1 0 

■:  004 

COMPLEX  F18.F1C.U2«U3«U4,U5.CMPLX 

lj  005 

DIMENSION  PH (43) *ZH(43)*XT( 10)«AT( 

X*  006 

DIMENSION  Z2( 10) .8J(50) 

007 

MP«NP-l 

!--■  008 

MT=  MP—  1 

009 

NxMT+MP 

010 

N2=2«N 

I  Oil 

00  11  K= 1 *NF 

"  012 

X=THR(KJ 

'  013 

CS(K)sCOSCX) 

014 

SMK)=SIN(X) 

■  015 

11 

CONTINUE 

-  016 

U=  (  0  .  •  1  •  ) 

1  017 

Ul=3. 141593*U**M1 

T  016 

M3=M1 ♦ 1 

oio 

M4»M243 

V  020 

IF(Ml.EO.O)  M3=2 

021 

M5=Ml*f  2 

:  022 

M 6s M2 42 

S  °23 

DO  12  I P— 1 • MP 

024 

K2=  IP 

. :  025 

I»IP*1 

•  026 

DR=.5*(RH( I ) — RH( IP) ) 

>-/  027 

DZ=.5*<  2H(  I  )-ZH(  IP)  ) 

Jv  028 

DI=SQPT ( DP*DRtDZ*DZ ) 

l029 

IF(IP.EQ.l)  RD=1 ./(01*Dl ) 

‘  030 

R3=«5*(RH(I )4RH{ IP)) 

031 

DR3xOR*R3*RD 

032 

RI=.5*R3 

033 

08=-02*R14RD 

^034 

D1R=DI*RD 

1  035 

D6=R3*DIR 

.  036 

21=.5*( 2H€ I )4ZH( IP) ) 

037 

0R=  »5*0R 

■;  038 

02=CR/Rl 

;  039 

DC  13  L=l.NT 

^040 

R2(L)=R1+0R*XT(L) 

*?041 

Z2(L)=Zl+0Z*XT(L) 

042 

13 

CONTINUE 

;*•  043 

DO  14  K=1,NF 

044 

CC=CSCK) 

045 

SS=SN(K) 

l'^046 

03=  OR  *CC 

j!o47 

04=-0Z*SS 

■  048 

07=-D6*CC 

’  049 

D9=-De*CC 

050 

D5=0R34SS 

.  051 

DO  23  M=M3  »M4 

a052 

FA<M>=0. 

rr-*053 

F8(NM0. 

;  054 

FC(M)*0* 

055 

23 

CONTINUE 

’*'056 

DO  15  L* 1 »NT 

‘-.057 

X=SS*R2CL) 

*056 

lF(X.GT..5E-7)  GG  TO  19 

—059 

00  20  M*M3 • M4 

>060 

8J(M)=0. 
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061 

20 

CCNT1NUE 

062 

BJf 2)-t • 

063 

S-=  1  . 

064 

GO  TO  18 

1  065 

19 

M=2.8*X*14.-2./X 

|  066 

I F( X  «LT • *5 )  M— 1 1 • 8+  AL OG 1 0 ( X ) 

067 

IF {  M»LT «M4  )  M=M4 

066 

BJt *1=0. 

069 

JH=M-1 

070 

BJI JH)=1. 

071 

DO,  16  J«4  •  M 

R  072 

J2— JM 

073 

JM= JM-l 

074 

Jt=JM~l 

L-:  075 

BJI JM|=J1/X*BJ I J21-8JI JMF2I 

M  076 

16 

CONTINUE 

i  077 

S-0. 

1  078 

IFfH.LE.4l  GO  TO  24 

079 

OC  17  J=4  • M  #  2 

C  060 

S=S*BJ| J1 

081 

17 

CONTINUE 

082 

24 

S=BJI2)«-2.*S 

■  083 

18 

AR G-Z 2 1  L  > *CC 

1  084 

UA=AT(L )/S*CMPLX{COS( ARG) 

• SIN(ARG) ) 

::  085 

UE=XT|L )*UA 

086 

UC=  XT (L )*UB 

087 

00  25  M=M3.M4 

088 

FA| H|=B JIM) *UA*FA(MJ 

121 

UC=I1./IM-2J>*UA 

M  089 

FB|Ml=BJ|Ml*UB*FBIM) 

1 22 

U5=C8*UC 

P  090 

FCIM1=BJ|M1  *UC*FC  1  M  | 

123 

U2*L5*l F2AfD2*F2B» 

.  091 

25 

CONTINUE 

124 

U3=U5*|F2B*02*F2CI 

y.  092 

15 

CONTINUE 

125 

U5=054UC 

;•  093 

IFIM1.NE.0)  GO  TC  26 

126 

UC=094UC*U 

/  094 

FA| ll=-FA|31 

127 

U4=UC*I  FIA4024F1B)*U5*<FACM)402*FB«M1> 

|  095 

FBI  1  J=-FB13» 

128 

U5=UC* 1 F1B*02»F1C 1*U5* 1 FBI  *  1 ♦02*FCI*O  1 

P  096 

FCI 1 )=-FCl3 1 

129 

R |K2*MT  1  =D7*UA*I  F2A*02*F281 

-  097 

26 

U  A=  U 1 

1  70 

R  |K5«-MT  )-=  D6  *UB*  1  F  1  A*-D2*F1 81 

098 

00  27  M=MS»M6 

131 

29 

1FI IP.EQ.l 1  GO  TC  21 

:■  099 

M7=H-l 

132 

R(K1>=R1K1 > ♦U2—U3 

100 

M8=MM 

1  33 

RIK41 =R IK41FU4-U5 

%  101 

F1A=FAI M8I-FAIM7) 

134 

IFI1P.EQ.NP)  GO  TO  22 

~  102 

F  t  B=F  B I  M8  1  —  F6 1  M7  1 

135 

21 

R<K21=U2+U3 

103 

UB=U*UA 

136 

RlK5l=U4+U5 

104 

K1=K2-1 

137 

22 

K2-K2+N2 

S  105 

K4=K1  *N 

138 

UA=UB 

106 

K5=K2+N 

139 

27 

CONTINUE 

^  107 

IFIM.NE.2)  GO  TO  28 

1  40 

1  4 

CONTINUE 

<r  108 

UC-034UB 

141 

12 

CONTINUE 

109 

U5=04*U  A 

142 

RETLRN 

V  110 

U2=UC*FIA+U5*FA( Ml 

l  43 

ENO 

:::  in 

U3=UC*F lfl«-U5*FBt  N 1 

C  112 

U4=0  • 

H  113 

U5*0« 

r  ti4 

R|K2*NT1=0. 

-  115 

RIK5FNT  1=0  1  R*U0*  IF l A* 02*F 1 B 1 

I  116 

GO  TO  29 

117 

28 

F2A=FAI N8 ) *F A 1 M7 1 

:  118 

F28=F8< *81 ♦FBIM7 ) 

119 

F2C=FCI M8JFFCIM71 

"  120 

F1C=FCIW81-FC(M7  > 

ii- 


» 


I 
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V.  THE  SUBROUTINES  DECOMP  AND  SOLVE 

The  subroutines  DECOMP (N,  IPS,  UL)  and  SOLVE(N,  IPS,  UL,  B,  X)  solve 
a  system  of  N  linear  equations  in  N  unknowns.  The  input  to  DECOMP  consists 
of  N  and  the  N  by  N  matrix  of  coefficients  on  the  left-hand  side  of  the 
matrix  equation  stored  by  columns  in  UL.  The  output  from  DECOMP  is  IPS 
and  UL.  This  output  is  fed  into  SOLVE.  The  rest  of  the  input  to  SOLVE 
consists  of  N  and  the  column  of  coefficients  on  the  right-hand  side  of  the 
matrix  equation  stored  in  B.  SOLVE  puts  the  solution  to  the  matrix  equa¬ 
tion  in  X. 

Minimum  allocations  are  given  by 
COMPLEX  UL (N*N) 

DIMENSION  SCL(N),  IPS(N) 
in  DECOMP  and  by 

COMPLEX  UL (N*N) ,  B(N) ,  X(N) 

DIMENSION  IPS (N) 

in  SOLVE. 

More  detail  concerning  DFCOMP  and  SOLVE  is  on  pages  46-49  of  [7]. 
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LISTING  OF  THE  SUBROUTINES  OECOMP  AND  SOLVE 
SUBROUTINE  OECOMP ( N  *  I  PS • UL ) 

COMPLEX  UL( 1600)  .PI VOT.EM 
DIMENSION  SCLI40 ) • I PS (40) 

DC  5  1=1  .N 
IPS!  1 1*1 
RN=0  • 

J1=I 

DO  2  J=1.N 

ULM=ABS ( REAL (UL ( Jill) 4ABS( A I MAG ( UL(  Jl 11) 

J  1=  J1  *N 

IF(RN-ULM)  1.2.2 

1  RN=ULM 

2  CONTINUE 
SCL (  I  )=  1 «/RN 

S  CONTINUE 
NMI=N-1 
K  2=0 

OO  17  K=1.NM1 
8IG=0. 

00  11  I =K  »N 
IP=IPSC II 
IPK=(P4K2 

S I ZE-  ( ABS ( REAL (UL( IPK )  ) ) +A8S ( A  I MAG( UL ( IPK ) ) >  > *SCL( IP) 
IFISI ZE-BIG)  11.11.10 

10  BIG=S IZE 
IPV=I 

11  CONTINUE 
IFIIPV-K)  14.1S.14 

14  J= IPS (  K  ) 

IPS(K)=IPS(  IPV) 

IPS(IPV)=J 

15  KPP=I PS( K )  +K2 
PIVOT=UL(KPP) 

KF1=K+1 

00  16  I=KP1.N 
KP=KPP 
IP=IPS(  IM-K2 
EM=-UL( IPl/PIVOT 


18  UL(IPI=-EM 

OO  16  J=KP1.N  061 

IP=IP4N  062 

KP=KP*N  063 

ULl IP)=UL«IP)*EM*UL(KP>  364 

16  CONTINUE  065 

K2=K24N  066 

17  CONTINUE  067 

RETURN  °68 

ENO  069 

SUBROUTINE  SOLVE ( N . I  PS . UL . 0 . X )  070 

COMPLEX  UL< 1600 ) .8(40 )• X< 40 ).  SUM  071 
OIMENSION  IPS( 40 )  072 

NP1*N41  073 

IP=IPS(1)  074 

X(  1  )=B( IP)  075 

DO  2  1=2  .N  076 

IP«IPS(t)  077 

IP8=  IP  078 

IMl»l-l  076 

SUM=0.  O8® 


00  1  J= 1 • IM 1 
SUM=SUM4UL( IP)*X(J) 

1  IP=  IP4N 

2  X(  1  )  =  B(  IPB )— SUM 
K2=N*(N-l ) 

IP= iPSt  N)4K2 
X ( N )  =  X ( N)/UL( IP) 

00  4  IBACK=2.N 
1  =  NP 1 - 1  BACK 
K2=K2-N 
IPI=IPS( I ) 4K2 
I P 1=1 ♦ 1 
SUM  =  0  • 

IP=IPI 

DC  2  J= IP  1 . N 
IP=IP+N 

3  SUM— SUM*UL (IP) *X ( J ) 

4  X( I )= ( X ( I ) -SUM ) /UL( IPI ) 
RETURN 

ENO 
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VI.  THE  MAIN  PROGRAM  FOR  THE  NEW  E-FIELD  SOLUTION 

The  main  program  for  the  new  E-field  solution  uses  the  subroutines 
ZMAT,  PLANE,  DECOMP,  and  SOLVE  to  calculate  the  surface  density  of* electric 
current  and  electric  charge  induced  on  S  by  an  axially  incident  plane 
wave.  For  the  surface  S  of  revolution  illuminated  by  the  9-polarized  plane 
wave  (60),  [1,  Eq.  (46)]  specializes  to 


.  v  Tm0  Tm  c  Teo  ,  Te 

J=  )  ()  I.J.+  )  I.kp.J.) 

—  L  nj  — nj  ,L,  nj  'j  — nj 

n=-°°  j  =  l  J  J  3=1  J  J 


(117) 


tn  6 

where  the  expansion  functions  J  .  and  kp.J  .  are  given  by  [1,  Eqs.  (91),  (92), 

-~nj 

(100),  and  (101)].  According  to  (90),  the  coefficients  I™®  and  I®®  in  (117) 

fr-Q 

are  the  elements  of  the  column  vector  I  that  satisfies 

n 


Z  I6  =  v9  ,  n  =  0,  ±1,  ±2,, 
n  n  n 


(118) 


where 


.0  J 


(119) 


,  n  =  ±1,  ±2, 


(120) 


(121) 


^  =  B 

n 


n  =  ±1,  ±2, 


(122) 


ft 
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It  is  evident  from  (15),  (2),  and  (13)  that 


z  =  8 

-n 


^mm  _2me 
n  n 


L„em  ee 

■Z  Z 

n  n  .. 


Equations  (87)  and  (122)  imply  that 


,  n  =  1,2,. 


(123) 


v9  =  8 

-n 


,  n  1,2,... 


(124) 


In  view  of  (123)  and  (124),  comparison  of  (118)  with  (118)  with  n  replaced 
by  -n  reveals  that  ’ 


,  n  1,2,... 


(125) 


It  is  assumed  that  the  angle  0t  of  incidence  of  the  plane  wave  is  either 
zero  or  tt  radians.  Consequently,  (98),  (105),  and  (113)  predict  that 


V  =  0,  n  i*  n 
n 


Substitution  of  (126)  into  (118)  gives 


T  =  0,  n  +  ±1 
n 


(126) 


(127) 


Equations  (125)  and  (127)  reduce  (117)  to 


where  I™?  is  the  jth  element  of  T™9  and  I^9  is  the  jth  element  of  I^9, 

y  1 1 1 Q  y 

According  to  (120),  the  combination  of  1^  and  1^  is  the  solution 
t9  of  (118)  for  n  =  1. 


(128) 
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The  expansion  functions  in  (128)  are  given  by  [1,  Eqs.  (100)  and 

(101)]  where  Jt.  and  J^.  are  given  by  [1,  Eqs.  (82)  and  (83)].  As  a 
— nj  -nj 

result,  (128)  reduces  to 


m  p 

J  =  u  J  k  cos  d)  +  u,  (J ,  +  J,)  k  sin  <p 
—  — t  t  — q>  <p  <J> 


(129) 


where 


P-2  T  (t) 

j  =  2  y  — )  i“? 

t  .£x  kp  lj 


(130) 


=  2  l  ( 


P-2  P.+1(t)  P.(t) 

'kA .  .  kA.  ’  lj 

J-l  .1+1  J 


(131) 


^  v*> 


(132) 


In  (129),  J™  is  due  to  the  magnetostatic  expansion  functions  J^\  and 

6  6 
in  (128),  and  is  due  to  the  electrostatic  expansion  functions  J  and 

J6^  in  (128).  Expression  (131)  can  be  rewritten  as 


where 


P-1  P.(t)  a 

jm  =  2  y  -j _ dm9  -  im0) 

<t»  jii  kA.  Ul,j-1  l]j) 


^-0 


(133) 


lmQ  =  0 
l.P-l 


(134) 


The  electric  charge  is  given  by  [1,  Eq.  (1)]  where  Vg*  and  J 
are  given  by  [1,  Eq.  (8-3)]  and  (129),  respectively.  Thus, 


1.  ’  j  5F  <°  \k  ♦>  +  H  (<JJ  +  J*>  k  sln  ♦» 


Thanks  to  expressions  ( 1 30)  —  ( 132 )  for  .Jt,  .J™,  and  J?,  (135)  becomes 


(135) 
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The  main  program  for  the  new  E-field  solution  is  listed  at  the  end 
of  this  section.  In  this  main  program,  input  data  are  read  from  punched 
cards  according  to 

READ (1 , 15)  NT,  NPHI 

15  FORMAT (2 13) 

READ (1 , 10) (XT (K) ,  K=l,  NT) 

READ (1,10) (AT (K) ,  K=l,  NT) 

10  FORMAT (5E14. 7) 

READ ( 1 , 10) (X (K) ,  K=1 ,  NPHI) 

READ(1,10)(A(K),  K=1 ,  NPHI) 

READ (1 ,  16)  NP,  BK,  THR(l) 

16  FORMAT (13,  2F.14.7) 

READ (1,18) (RH(I) ,  1=1,  NP) 

READ ( 1 ,  18)  (7,H(I) ,  1=1,  NP) 

18  FORMAT (10F8. 4) 

NT,  NPHI,  XT,  AT,  X,  and  A  are  Gaussian  quadrature  data  that  are  fed  into 

the  subroutine  7.MAT.  NT  is  n^  in  [2,  Eqs.  (62 )  —  (63)  ] .  The  variable  n^ 

that  appears  in  [2,  Eqs.  (132)— (133) ]  was  chosen  to  be  equal  to  n  .  Hence, 

<\> 

NT  can  also  be  viewed  as  n^.  XT  contains  the  n^  abscissas  x^,  in 

(nt> 

[2,  Eqs.  (62)— (63) ] ,  and  AT  contains  the  nt  weights  A^,  in  [2,  Eqs. 

(62)— (63) ] .  NPHI  is  in  [2,  Eqs.  (64 ) — (66 ) ] .  X  contains  the  n^  abscissas 

(n  )  (n  ) 

Xp  ^  of  [2,  F.q.  (70)],  and  A  contains  the  n^  weights  A^  p  in  [2,  Eqs. 

(64)-(66)].  NP  is  the  number  P  of  data  points  on  the  generating  curve  of  S. 
P  appears  in  [1,  F.qs.  (100)  and  (101)  j.  BK  is  the  propagation  constant  k 
in  1/meters,  k  appears  in  (140)  and  (141).  THR(l)  is  the  angle  of  inci¬ 
dence  6  of  (62)— (63)  that  is  fed  into  the  subroutine  PLANE.  THR(l)  is  in 


radians.  THR(l)  should  be  either  zero  or  it.  If  THR(l)  is  zero,  the 
incident  electric  field  is  given  by  (140).  If  THR(l)  is  tt,  the  incident 
electric  field  is  given  by  (141).  RH  and  ZH  contain  the  coordinates 
of  the  data  points  on  the  generating  curve  of  S. 

RH(J)  =  p(t“),  .T  =  1,2, . .  .NP  (142) 

ZH(J)  =  z(t“),  -T  =  1,2, ...NP  (143) 

Here,  z  is  the  coordinate  along  the  axis  about  which  the  generating  curve 
of  S  is  rotated,  and  p  is  the  distance  from  this  axis.  Also,  (t  ^)  denotes 
evaluation  at  the  Jth  data  point.  Tn  (142)-(143),  p(tj)  and  z(t^)  are  in 
meters.  The  main  program  for  the  new  E-field  solution  uses  the  subprograms 
ZMAT ,  BLOG,  PLANE,  DECOMP,  and  SOLVE  of  Sections  II,  III,  IV,  and  V. 

Minimum  allocations  in  this  main  program  are  given  by 

COMPLEX  Z (N*N) ,  R(2*N) ,  B(N) ,  C(N),  CE(NP-l), 

CQ(NP-l) 

DIMENSION  RH(NP) ,  ZH(NP),  X(NPHI),  A(NPHI), 

XT (NT),  AT (NT),  IPS(N) 


where 


N  =  2*NP-3 


(144) 


With  regard  to  (129)  and  (136),  the  main  program  prints  under  the 
heading  "T  COMPONENT  OF  ELECTRIC  CURRENT,"  +  J*  under  the  heading  "PHI 

COMPONENT  OF  ELECTRIC  CURRENT,"  under  the  heading  "MAGNETOSTATIC  PART  OF 

$ 

JP,"  J®  under  the  heading  "ELECTROSTATIC  PART  OF  JP,"  and  q  under  the  head- 

<J> 

ing  "ELECTRIC  CHARGE."  A  real  part,  an  imaginary  part,  and  a  magnitude 
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*  =  Vl 


(145) 


tn  6  m  6 

The  jth  line  of  printed  output  for  each  of  +  J^,  J^,  J^,  and  q 

is  the  value  at  the  center  of  the  domain  of  (t) .  This  center  occurs 


‘  '  2  (t3  +  Vl* 


(146) 


The  sample  input  and  output  data  are  for  the  conducting  circular 

disk  of  radius  0.002  wavelengths  illuminated  by  the  incident  electric 

field  (140).  The  disk  is  placed  in  the  xy  plane  and  is  centered  at  the 

origin.  The  factor  k  in  (129)  and  (136)  compensates  for  the  factor  k  in 

m  g 

(140)  so  that  J ^ +  J^,  and  q  coincide  with  the  variables  J^,  J^,  and 
q  used  in  [1,  Eqs.  (108)  and  (109)].  The  printed  output  for  j  J  ^  j  is 
plotted  by  means  of  the  interior  x's  in  [1,  Fig.  1].  The  printed  output 
for  |j™  +  J^|  is  plotted  with  the  x's  in  fl.  Fig.  2].  The  printed  output 
for  | q |  is  plotted  with  the  x's  in  [1,  Fig.  3]. 

DO  loop  28  multiplies  RH(J)  and  ZH(J)  by  the  propagation  con¬ 
stant  k.  With  regard  to  (139),  line  41  puts  7.^  in  Z.  Line  44  prints 
out  the  first  column  of  Z^.  Line  46  puts  of  (122)  in  R(l)  through 
R(N)  where  N  is  given  by  (144).  Line  46  also  puts  the  column  vectors 

and  BVe^  of  (89)  in  R(N+1)  through  R(2*N),  but  these  column  vectors 

are  not  useci  in  the  main  program.  DO  loop  22  and  line  52  take  advantage 

->0  ->-0 

of  (124)  to  store  in  B ( 1 )  through  B(N).  is  needed  because  it 
appears  in  (139).  Lines  55  and  56  put  the  solution  1^  of  (139)  in  C(l) 
through  C (N) . 

DO  loop  24  prints  out  J  of  (130).  Inside  DO  loop  27,  line 
85  puts  of  (132)  in  CE(.T),  line  86  puts  q  of  (137)  in  CQ(.I),  and 


line  90  prints  out  4-  and  of  (133).  As  intermediate  steps, 
line  79  puts  2/(kA^)  in  C4,  lines  80-83  put  J™  of  (133)  in  Cl,  and 

m  0 

line  87  puts  in  C3.  Inside  DO  loop  32,  line  103  prints  out 


-ISTING  OF  THE  MAIN  PROGRAM  FOR  THE  NEW  E-FIELD  SOLUTION 
THE  SUBPROGRAMS  ZMAT,  BLOG*  PLANE*  DECOMP*  AND  SOLVE  ARE  NEEDED 
JOB  ( XX  XX • XXXX • 1 • 1 1  *  *  MAUT  Z* JCE  *  «  REG IQN=  20 OK 
C  WATFIV 

rsiN  oo  ♦ 

MAUTZ  *T I  ME -5  *PAGES=60 

COMPLEX  ZI 1600I.RI240 )«Bt  40) *C(40)*U*CI *C3*CE< 20) *CO(20 » 

>  I  MENS  I  ON  THRO)  *RH(  4  3)  •  ZHI  431  •  X(  48) • A< 48) . XT! 10) • AT! 10) • IPSI40) 
READ! 1*15)  NT«NPHt 
=CRMATI 213) 

•RITE ( 3*30 )  NT.NPHI 

FORMAT ( •  NT  NPH I •/ 1 X  *  13* 15 ) 

READ! 1* lO)IXT(K) ,K=1*NT) 

RE ADI  1 • 10)(AT(K) *K=1,NT) 
r CRM AT  I 5E1 4 • 7 ) 

•  RITEI3*U)tXTIK).K=l .NT) 

•RITEI3*12)t AT ( K )  •  K  =  1 *NT ) 

=ORMAT I •  XT*/I IX *SE14.7) ) 
rGRMAT I •  AT*/I IX  .SE14.7) ) 

READ! 1*  10) ( X( K) , K=1.NPHI ) 

READ!  1  •  10)1 A|K) *K=1 *NPHI  ) 

•RITE <3. 13) (XI K) .K=1*NPHI) 

•  RITE 13*14)1  AIK) *K=l*KPHt ) 

-ORMATI*  X • ✓< 1 X*  5E 1 4*  'I ) 

=  0RMAT I •  A* / 1  1 X  *  SE 14*7) ) 

READ!  1*16)  NP#tlK  •  THRI  1 ) 

-CRMATI I3.2E14.7  ) 

•  R  ITE  13*17)  NP.BK.THRU) 

:0RMAT I •  NP* *6X • *  BK *  * 12X*  * THR* /I X* I  3 • 2€ 14*7) 

IE ADI  1 • 18) IRHI I ) *1  =  1 *NP) 

READ! 1*18)1 ZH I  I) *1  =  1* NP ) 

'ORMATI 10F8.4) 

(RITE  13. 19) IRHI I ). 1  =  1 *NP) 

(RITE  I  3* 20) I ZHI I  ).I  =  1«NP) 

CRMATI*  RH*/I1X*10F8.4)) 
rCRMAT I •  ZH*/I IX* 10F8.4) ) 

>0  28  J=1*NP 
1HI J)=BK*RHI J) 

!Ht  J)=BK4ZHt J) 

:cntinue 

:ALL  ZMATI 1  *  1 *NP.NPHI *NT*RH.ZH.X* A*  XT* AT.Z) 

(T=NP— 2 
1=2  *MT  ♦  1 

i  RITE  13.29) I Zt J)  *J=1.N) 

CRMATI*  Z*/I1X.6E1 1*4)1 

ALL  PLANE!  1*1.1  .NP . NT , RH, ZH. XT, AT* THR* R » 

10  22  J=1.MT 
I J)=RI J) 
lsJtMT 
I J1 )=— RI J1 I 
ONTINUE 
I N )=— R IN) 

RITE  13 .23) (BIJ) • J=1.N» 

ORMATI*  B*/tlX.6E11.4)> 

ALL  DECOMPIN.IPS.Z) 

ALL  SOLVE IN*IPS*Z*B»C) 

RITE  13*37) 

ORMATI *0  T  COMPONENT  OF  ELECTRIC  CURRENT*) 

RITEI3.21 ) 

ORMATI*  REAL  JT 


IMAG  JT 


MAG  JT*) 
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MAG  JP*.6X. 'REAL 


IMAG  JP 


>1  00  24  J= 1 • MT 

♦2  C1=2./RH<  J41  I4CC  Jl 

>3  C2=CABSIC1» 

*4  WR  I TE ( 3 • 25 )  C1.C2 

Is  25  FORMAT ( 1X.3E11.4) 

>6  24  CONTINUE 

>7  U=<0..2.» 

.8  HR  I  TE  (  3*26  ) 

>9  26  FORMAT ( • OPH I  COMPONENT  OF  ELECTRIC  CURRENT* • 7X. • MAGNETOSTATIC  PART 

fO  1  OF  JP*  > 

ll  WRITE ( 3*351 

2  35  FORMAT { •  REAL  JP  IMAG  JP  MAG  JP*.6X.*REAL  JP  IMAG  JP 

1  MAG  JP*I 
MPsKP-1 
DO  27  J*  I  #MP 
J  F— J4 1 

D 1=RH< JPl—RHf  J ) 

D2=ZHC JP )  — ZH<  J ) 

C 4= 2. /SORT! D14D1 402402) 

0=0. 

IF(J.NE.l)  Cl=C14CIJ-l» 

IFfJ.NE.MP)  C1=C1-C(J) 

C t=C4  *C 1 
C3=C4 J4MT) 

CE( J|=U*C3 

CQ4 J)=-4./4RH4 JP14RH4 J»»*C3 
C3=C1*CE4 Jl 
C4=CA8S4C3» 

C2=CAES(C1 I 

WR  I  TE 13*33)  C3.C4.C1.C2 
33  FCRMATI 1X.3EI1.4.2X.3E1 1.4) 

27  CONTINUE 

WR I TE43.31 ) 

31  FORMAT I*0*.5X.' ELECTROSTATIC  PART  OF  JP * . I 5X . • ELECTR IC  CHARGE* I 
WR I TE (3 .36 1 

36  FORMAT! •  REAL  JP  IMAG  JP  MAG  JP*.6X.*REAL  O  IMAG  O 

1  MAG  0* ) 

OO  32  JxI.MP 
C1=CE<  J» 

C2=CABS4C1» 

C3=C0( J) 

C4=CABS4C3 ) 

WRITE43.331  CI.C2.C3.C4 

32  CONTINUE 
STOP 
ENO 

SOATA 
2  20 

•0.5773503E400  0.5773503E400 
O.IOOOOOOE+O 1  0.  100000CE401 

0.9931286E400-0.9639719E400-0.9122344E400-0.8391 170E400-0.746331 9E4  00 
4-0.636053 7E400-0.510867CE*  00- 0. 3 73706 1 E400-0. 2277*S9E*00-0 .T652652E-01 
0. 7652652E-0 t  0. 2277859E400  0. 3737061 E400  0. 51 08670E400  0 • 636053 7E400 
0. 746331 9E 400  0 .8391 1 70E400  0 . 9 122 344E 40 0  0.963971 9E400  0. 9931 286E400 
0. 1 761401E-01  0 .4060 1 4 3E  — 0  I  0.6267205E-0 1  0. 8327674E-0I  0 . 1 0 1 9 30  IE 400 
0. 1 181945E 400  0. 13168eeE400  0.  1420961E400  0 • 1 49 1 730E 400  0 • 1 527534E*00 
0. 1527534E400  0. 1491 730E400  0 . 1 42096 1 E 400  0. 1 316886F400  0. 1  1 81 945E400 
4 0. 1C19301E400  0 • 83276 7 4 E— 01  0. 62672 05E-0 1  0 • 4060 1 43E-0 1  0 • 1 76 140 1E-0 1 
16  0. 83  7  75  COE— 03  0. 00 C0000E400 

0.0000  1.0000  2.0000  3.0000  4.0000  5.0000  6.0000  7.0000  8.0000 


4.0000 


>0000 


>0000 


.0000  8.0000 


.0000 
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% 


10.0000 
0.0000 
0.0000 
•  STOP 
* 

✓ 


11.0000 

0.0000 

0.0000 


12.0000 

o.oooc 

0.0000 


13.0000 

0.0000 

0.0000 


14.0000 

0.0000 

0.0000 


15.0000 

0.0000 

0.0000 


0.0000  0.0000  0.0000  0.0000  \ 


»R 1 NTEO  OUTPUT 
NT  NPHI 
2  20 

r 

®L0.  5773503E  400  0. 


5773503E+00 


:*T 

.0.10000006+01  0.10003006+01 

r-’x 


>-0.9931 28CE+00— 0. 
BO.  636053  7E  +  00-0. 
0.  7652652E-01  0. 
0. 7463319E+00  0. 

0. 1761401E-01 
: 0. 1181945E+00 

BO.  1527534E+00 
0. 1019301E+00 
NP  BK 

^  16  0.83 7 75  COE— 03 

S  0.0000  1.0000 

10.0000  11.0000 

M 

0.0000  0.0000 
0.0000  0.0000 


0. 

0. 

0. 

0. 


96397 19E +00-0. 9 122344E+0 0-0. 8391 1 70E  +  00- 
5 1  086  70E+00— 0 .37 37061E+00-3.2277 859E+00- 
22778S9E+00  0 . 373706 1  E  +  00  0.  51  086  70E4-00 
8391 1 7064-00  0.91 223 44 E 4 00  0 .96397 l 9E+00 

4060143E-01  0 • 62672 08E—0 1  0 . 83276 75E-0 1 
13 168C6E+  00  0. 1420961E+00  0 . 1 49 1 730E+00 
1491720E+00  0.14209616+00  0.13168866+00 
8327675E-01  0 .6267208E-0 1  0. 4060 1 43E-0  1 
THR 

0.00C0000E+00 

2.0000  3.0000  4.0000  5.0000 

12.0000  13.0000  14.0000  15.0000 


0. 746331 9E  4-00 
0.76S2652E-01 
0. 636053 7E400 
0.9931  2866+00 

0. 1 019301E400 
0 . 1527534E400 
0.1  181945E+00 
0. 1 761 401E— 01 


rl 


6.0000  7.0000  8.0000  9.0000 


0.0000 

0.0000 


0.0000 

0.0000 


0.0000 

0.0000 


0.0000 

0.0000 


0.0000  0.0000  0.0000  0.0000 


_ f 0. 4101 E—  04 
|jo.3052E-04- 
^0.00006+00 
0.3052E-  C4" 
0.  12216-03- 
0.7732E4-01 
0.  7767E+00 
0.2670E4-00 
^0.13466+00 
0.8102E— 01 

0.00006+00 
0. OOOOE  400 
®O.OOOOE+O0 
0. OOOOE4QO 
O.OOOOE+OO 
O.OOOOE+OO 
0. OCOCE+OO 

.o.ooqoe+oo 
•b.oooo E+OO 
0.00006+00 


0.2046E+05- 
-0.529EE+03- 
-0.8246E+02 
-0.272IE+02- 
-0.1220 E  +  02- 
0.2608E-07 
0 .40S6E—  07 
0.6333E— C7 
0.694  IE— 07 
0.1192E-06 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

0*1  885E+02- 

0  .5655E+02- 

0.94256+02- 

0.131 9E  +  03- 

0.1696E+03- 


-0.37196—04' 
-0.4578E— 04- 
0.61 04E-04 
-0.15266-04 
-0.91  556-04- 
0.31 26E+01 
0.5074E+00 
0. 20676  +  90 
0.  U20E  +  00 
0.70166-01 

O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
■C. OOCOE+OO 
■O.OOOOE+OO 
■O.OOOOE+OO 
■O.OOOOE+OO 
•O.OOOOE  +  OO 


-0.32266+04 
-0  •  2448E+0  3- 
-0.5424E  +  02- 
-0.20306+02- 
-0  •  9734E  +  0  1- 
0.4098E-07 
0.7078E— 07 
0.00006+00 
0.1006E-06 
0.17886-06 

O.OOOOE+OO 
O.OOOOE+CO 
O.OOOOE+OO 
0 .00006+00 
O.OOOOE+OO- 
0.31 426+02- 
0.691 2E+C2- 
0.  10686  +  03- 
0. 14456+03- 
0. 1822E  +  03 


0.15266-04 
0.457BE-04 
0. 1526E-03 
0.1526E-03 
0.36506+01 
0.13536+01 
0.3584E+00 
0 .16496+00 
0.9462E-0 1 


0.00006+00 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
0. OOCCE+OO 
O.OOOOE+OO 
O.OOOOE+OO 


-0.1689E+04 
-0.1  347E+03 
-0.3764E+02 
-0. 1558E+02 
0.861 5E -08 
0.1490E-07 
0 .81966-07 
0.12296-06 
0.2384E— 06 


O.OOOOE+OO 
0*00006+00 
O.OOOOE+OO 
0. OOOOE +00 
0 .62836+01 
0 .4398E+02 
0.8168E+02 
0.1 194E+03 
0.15716+03 


t  ' 


f  COMPONENT  OF  ELECTRIC  CURRENT 
REAL  JT  IMAG  JT  MAG  JT 

j. 75136-08  0.21  16E-01  0.2116E-01 
5.71336-08  0.2095E-01  0.20956-01 
>•82266-08  0 .20696-01  0.2069E-01 

IT  *!■'.  'tm.  .  i  -  ml  fcLai 


jr  T» 


0.8463E-08  0.2034E-01  0.2034E-01 
0« 8348E— 08  0 • 1988E-0 1  0.1988E-01 
0.7887E-08  0.1931E-01  0.1931E-OI 
0.7588E-08  0.1862E-01  0.1862E-0I 
0.7201E-08  0*1 778E-01  0.I778E-0I 
0.  66  82  E- 08  0.1679E-0I  0.1679E-01 
0.6318E-08  0.1560E-01  0.1560E-01 
0.5688E-08  0.1418E-0I  0.1418E-01 
0.4960E-C8  0.1244E-01  0.I244E-01 
0.4I43E-C8  0.10 23E— 0 1  0.1023E-01 
0.2965E-08  0.7317E-02  0.7317E-02 

PHI  COMPONENT  OF  ELECTRIC  CURRENT 
REAL  JP  IMAG  JP  MAG  JP 

-0.7453E-08-0.2t 19E-01  0.21 19E-01 
-0.6943E-08-0.21 06E-01  0.21 06E-OI 
— 0. 1 079E— 07— 0.21 07E— 01  0.2107E-01 

—  0.1 020E— 07— 0 .2 107£— 0 1  0.2107E-01 
— 0. 9356E— 08— 0 .21 10E— 01  0.211 0E-01 
-0. 7953E-08— 0.21 12E-01  0.21 12E-01 
—0.91 1 4 E— 08— 0.2 1 20E— 0 1  0.21 20E-01 
-0. 1 004E— 07— 0.2130E-01  0.2130E-01 
-O. 8  794 E— 08  — 0 .21 46E— 0 1  0.2148E-01 
-0.1231 E-07-0.21 79E-01  0.21 79E-01 
-0. 1114E-07-0.2232E— 01  0.2232E-01 
-0. 1310E— 07— 0.2326E-01  0.2326E-01 
— 0 • 1 6 98E— 07-0 .25 1 2E— 0 1  0.2S12E-01 
-0. I 809E— 07— 0.2828E— 01  0.2828E-01 

—  0. 4S86E—  07—0 .S684E— 01  0.5684E-01 

ELECTROSTATIC  PART  OF  JP 
REAL  JP  IMAG  JP  MAG  JP 

0. 5986E— 10 -O .34 37E- 04  0.3437E-04 
-0. 1 9 l OE—  09—0.3 1 46E— 03  0.3146E-C3 
— 0. 3801E— 09— 0.900  IE— 03  0.90C1E-03 
-0. I025E-08-0. 1798E-02  O  .1798E-02 
-0.1 467E-08-0.3044E-02  0.3044E-02 
-0. 2371 E-C8-0.4667E-02  0 .4667E-02 
-0.331 8E— 08— 0.6742E— 02  0.6742E-02 
—0.5552 E— 08— 0.93 5 3E— 02  0.9353E-02 
— 0. 6258E— 08— 0 . 1265E-01  0.1 265E-01 
— 0. 9266E— 08— 0 • 1685E-01  0 .1685E-0 1 
-0.11 76E—  0  7—0 • 2236E—  01  0.2236E-01 
-0. 161 5E- 07-0 .2993E-01  0.2993E-01 

-0. 2263E— 07— 0.4143E— 0 1  C.41 43E-01 
—0.301 6E— 07 -0.5884E-01  0 .5884E-01 
-0. 8765E-07-0. 1593E400  0.1593E+00 


MAGNETOSTATIC  PART  OF  JP 
REAL  JP  IMAG  JP  MAG  JP 

-0. 7513E-08-0.21 16E— 01  0.21 16E-0I 
— 0.6752E— 08— 0.2075E— 01  0.2075E-01 
-0. 1041E-07-0.201 7E-01  0.201 7E-01 

—  0.91 72E— 08-0 . 1 92  7E— 0 1  0.1927E-01 
— 0. 7889E-08-0. 1 806E— 01  0.1806E-01 
— 0.5582E-08— 0. 1645E-01  0.1645E-01 

—  0. 5795E— 0  8—0 . 1 446E— 0 1  0.1446E-01 
— 0. 4488E— 08— 0. 1 1 94E— 0 1  0.1 194E-01 
-0.2536E— Oe— 0.8832E— 02  0.883 2E -02 
— 0.3040E— 08— 0.4939E— 02  0.4939E-02 

0.61 28E— 09  0.421 7E— 04  0.4217E-04 
0.3051E-08  0.6673E-02  0.6673E-02 
0.5654E-08  0.1632E-01  0.1632E-01 
0. 1 207E-07  0. 30S6E— 01  0.3056E-01 
0.4179E-07  0.1024E400  0.1024E400 

ELECTRIC  CHARGE 
REAL  0  IMAG  Q  MAG  O 

0.8204E— 01  0.1429E-06  0.8204E-01 
0. 2 50 3E 40 0—0 • 1 520E— 06  0.2503E+00 
0.4298E4C0-0. 1815E-06  0.4298E+00 
0.61 33E400-0 . 3496E— 06  0.6133F4-00 
0.8074C+00— 0.3891E— 06  0.8074F+00 
0. 1013E+01-0.5147E-06  0.1013E+01 
0.1 238E40 1-0.6093E— 06  0.1238E  +  01 
0.14896*01 — 0.8836E— 06  0.1489E*01 
0.1 7766*01-0. 8789E— 06  0.1776E*01 
0.21 17E*0 1-0.1 164E-05  0.2U7E*01 
0 . 2 54 2 E 40 1-0.1 336E— 05  0.2542E+01 
0.31 07E  +  0 l -0 • 1 677E— 05  0.3107E+01 
0.3957E*01— 0.21 61 E— 05  0.3957E+01 
0.5203E40 1-0.2667E-05  0.5203E*01 
0.131 lE*02-0. 72156-05  0.1311E*02 


VII.  BOUWKAMP'S  POWER  SERIES  SOLUTION  FOR  A  CIRCULAR  DISK 


The  known  solutions  for  the  electric  current  and  electric  charge 
on  a  conducting  circular  disk  excited  by  an  axially  incident  plane  wave 
are  plotted  in  [1,  Figs.  1,2,  and  3].  The  known  solution  for  the  electric 
current  was  obtained  from  Bouwkamp's  formulas  [3,  Eqs.  (23),  (38),  (39), 
Table  I,  and  Table  II].  Actually,  Bouwkamp's  formulas  are  for  a  disk  of 
unit  radius.  They  had  to  be  changed  to  mksc  units  [4,  p.  1]  for  a  disk 
of  radius  a  meters  before  being  used.  The  electric  charge  was  obtained 
from  the  electric  current  via  the  equation  of  continuity.  The  foregoing 
electric  current  and  electric  charge  were  calculated  by  means  of  a  com¬ 
puter  program.  This  program  is  listed  at  the  end  of  this  section. 

In  the  present  report,  e~luJt  time  dependence  is  assumed.  However, 

Bouwkamp  uses  e  iU>t  time  dependence.  His  formulas  can  be  made  valid  for 

eJtut  time  dependence  by  replacing  i  by  -j  everywhere.  Henceforth  in  the 

present  report,  Bouwkamp's  formulas  will  be  taken  not  as  they  stand  but 

with  i  replaced  by  -j  everywhere. 

inc 

If  the  electric  field  E  given  by 


inc 


=  u  r)  e 


-jkz 


(147) 


is  incident  upon  a  conducting  circular  disk  of  radius  a  meters  lying  in 
the  xy  plane  and  centered  at  the  origin,  then  the  electric  field  integral 
equation  for  the  surface  density  J  of  electric  current  induced  on  the  disk 
is 

u  =  jfk  A (r )  +  ~  7 (V  •  A (r) ) ]  r  on  S  (148) 

x  k  r.  cin 


where 
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-  ta 


!  ff  i(k.r') 


(149) 


The  subscript  tan  denotes  the  component  in  the  xy  plane,  S  is  the  sur¬ 
face  of  the  disk,  and  ds'  is  the  differential  element  of  surface  area  at 
r' .  All  quantities  in  ( 147 ) — ( 149 )  are  in  mksc  units.  We  want  to  solve 
(148)  and  (149)  for  J. 


On  the  other  hand,  it  is  evident  from  [3,  Sec.  2]  that  Bouwkamp 


considers 


k  ■  J[k1A1<rl>  +  jq-  V’,  •  •  %  on  S1 


where 


=  C 


-jk, |r.-r{| 

1  ff  i(kr  ip  e 


r  -r' 

l£l 


(150) 


(151) 


He  solves  (150)  and  (151)  for  _I*  Here,  r^  and  r^  are  dimensionless 
radius  vectors,  k^  is  the  dimensionless  wave  number,  7^  is  the  V  operator 
with  respect  to  the  coordinates  of  r^,  ds^  is  the  differential  element  of 
surface  area  at  r^,  is  the  surface  of  a  disk  of  radius  unity,  and  c  is 
the  speed  of  light. 

Substitution  of 


-1  a 


(152) 


into  (150)  and  (151)  gives 


k^  =  ka 


(153) 


iL  =  j[ka  Ax(-)  +  £  V(V  *  X  °n  S 


1  a  tan 


(154) 


where 


kV  -  t 


-j  k | r-a  r  j| 

I(ka,  r!)  e 


(155) 
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In  (154),  V  operates  on  the  coordinates  of  r.  Introduction  of  a  new 
radius  vector  r '  defined  by 


r ' 


(156) 


reduces  (155)  to 


t  -  v  -jk|r-r' 
ff  I (ka,  — )  e  J  1 - 


A1  (-)  =  — 
— 1  a  ac 


r-r 


ds’ 


Equations  (154)  and  (157)  can  be  rewritten  as 


(157) 


u^  =  j[k  A(r)  +  i  V (V  •  A(r))]tan  ,  r  on  S 


(158) 


where 


4t r 


A(r)  = 


4tt 


I(ka,  — ) 


~ j k | r-r * 


r-r 


ds’ 


Equations  (158)  and  (159)  were  derived  from  (150)  and  (151)  for 
a  disk  of  unit  radius.  However,  (158)  and  (159)  will  coincide  with 
(148)  and  (149)  for  the  disk  of  radius  a  if 

J (k,  r)  -  I(ka,  j) 

C  3 


(159) 


(160) 


In  (160),  the  radius  vector  r_  on  the  surface  of  the  disk  of  radius  a 

is  specified  by  the  distance  p  from  the  center  of  the  disk  and  the  azimuthal 

angle  <J>  measured  from  the  positive  x  axis.  According  to  (160),  the  electric 

inc 

current  J  in  mksc  units  induced  on  the  disk  of  radius  a  by  E  of  (147) 
is  Bouwkamp's  solution  I_  for  the  disk  of  unit  radius  multiplied  by  4n/c 
with  Bouwkamp's  k  replaced  by  ka  and  with  Bouwkamp's  p  replaced  by  p/a. 

If 


J  =  u  J  +  u  J 
-  -xx  — y  y 


(161) 


then,  thanks  to  the  considerations  in  the  previous  paragraph  and  in  the 
second  paragraph  of  this  section,  Bouwkamp's  formulas  [3,  Eqs.  (23), 
(38),  (39),  Table  I  and  Table  IT]  yield 
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where 


where  the  first  six  A 


and  the  first  six  B  ' 


j  4(A  +  B  cos  2$) 

it  Vl  -  R2 

(162) 

j  4B  sin  2<J> 

y  Wl-R2 

(163) 

R  =  p/a 

(164) 

CO 

A  =  l  A  (-jka)n 
n=l 

(165) 

00 

B  =  l  B  (-jka)n 
n-1  n 

(166) 

's  are  given  by 
n 

-4  +  3R2 

A1  =  3 

(167a) 

a2  =  0 

(167b) 

56  -  40R2  +  5RA 

A3  =  90 

(167c) 

r2) 

(167d) 

-2656  +  2408R2  -  448R4  +  21R6 

A5  "  12600 

(167e) 

A6  =  6 JF  ("296  +  192r2  "  15r4) 

(167  f  ) 

s  are  given  by 

58 


B  -  — 
B1  "  3 


(168a) 


b2  =  o 

_  R2 (-8  +  R2) 
B3 - 30 


(168b) 


(168c) 


B  =  _  Ml 

B4  97t 

-  R2 (200  -  68R2  +  3R4) 
5  "  2520  . . 

_  _  2R2 (134  -  15R2) 


(168d) 


(168e) 


(168f ) 


As  an  alternative  to  (161)  -  (163),  J  can  be  written  as 


J  =  u  J  +  u.J, 
-  -DP  ^  <f) 


(169) 


where  and  are  the  unit  vectors  in  the  p  and  <j>  directions,  respec¬ 


tively,  and 


4  (A  4-  B)  cos 


it  VI  -  R 


(170) 


4(B  -  A)  sin  4> 

<J>  =  / - 2 

T7  VI  -  R 

It  is  evident  from  (165)  and  (166)  that 


(171) 


A  +  B  =  l  (A  +  B  ) (-jka)n 


(172) 


B-A=  )  (B  -A  )(-jka) 
L ,  n  n  J 
n=l 


(173) 


Equations  (167)  and  (168)  give 


ai  +  Bi  =  -  |  (1  -  r2) 


A2  +  B2  =  0 


A3  +  B3  =  ~  (!  -  R2) (7  -  R2) 


(174a) 


(174b) 


(174c) 
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C 

P 


+  c 


p2 


+  CP3  +  C 


p4  +  Cp5 


V 


V  +  Cd,2  +  C4)3  +  C4>4  +  C4>5 


where 


where 


Also, 


where 


j  4ka 
3 


S 

P 


r  -  j4(ka)3(7  -  R2) 
p2  “  45  Sp 

r  =  8^ka>A  c 
p3  9  it  bp 

C  _~j(ka)5(-  664  4-  188R2  -  9rS 
P 4  3150 

c  =  4(ka)6 (-  148  +  15R2) 

P5  675tt  Sp 


<,  =  4  /l  -  R2 

P  7T 

r  -  1 2ka (2  -  R2) 

V  3  S4> 

r  =  -i(ka)3(28  -  8R2  +  R4)  c 
<t>2  45 

r  =  8(ka)A  c 
<|>3  9tt  <jj 

C  =  -  J  Cka)5C-1328  +  704R2  -  54R4  +  3R6) 

4>4  6300  S<p 

_-4(ka)6(-148  -I-  29R2) 

<p5  675tt  Scj) 


IT 


4 


R 


2 


Substitution  of  (176)  and  (177)  into  (169)  gives 


(180) 

(181) 

(182a) 

(182b) 

(182c) 

(182d) 

(182e) 

(183) 

(184a) 

(184b) 

(184c) 

(184d) 

(184e) 

(185) 


J  =  u  C  cos  <j)  +  u.C.sin  (J) 
~  ~P  P  <p  * 


(186) 
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Thus,  the  electric  current  J  induced  on  the  disk  of  radius  a  by  the  inci- 

inc 

dent  electric  field  E  of  (147)  is  given  by  (186)  where  and  are 
given  by  (180)  and  (181),  respectively. 

The  electric  charge  associated  with  J  of  (186)  is  called  qg  and  is 
given  by  the  equation  of  continuity  [1,  Eq.  (1)].  In  view  of  [1,  Eq.  (B— 3) ] , 
substitution  of  (186)  into  [1,  Eq.  (1)]  gives 


-r~  [~  (pC  cos  (p)  +  (C  sin  <j>)] 
-.loop  dp  P  Y  d<j>  <|> 


(187) 


which  reduces  to 


C 

q 

q  =  cos  <h 
e  c  r 


(188) 


where 


C  =  (C  +  C  +  p 
q  kp  p  <p  dp 


(189) 


Thanks  to  (180)-(185),  (189)  becomes 


C  =C1+C0+C.+C.+Cc 
q  ql  q2  q3  q4  q5 


(190) 


where 


(191a) 


(ka) 2(4  -  K2) 


(191b) 


j 16 (ka) ~ 


(191c) 


(ka) ^ (88  -  44R2  +  3R4) 


(191d) 


jl6(ka)5(-26  +  5R2) 


(191e) 


where 


s 

q  ttVI  -  R2 


(192) 
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The  electric  current  .J  induced  on  the  disk  is  given  in  terms 

of  0^  and  by  (186).  The  electric  charge  qg  induced  on  the  disk 

is  given  in  terms  of  C  by  (188).  Now,  C  ,  C.,  and  C  can  be  cal- 

q  P  P  q 

culated  by  the  computer  program  listed  at  the  end  of  this  section. 
This  program  reads  input  data  from  a  punched  card  according  to 


READ (1,12)  N,  BK 
12  FORMAT (13 ,  E14.7) 

The  coefficients  C  ,  C,,  and  C  are  calculated  and  printed  out  at 

P  <J>  q 


p/a  =  0, 


1  2 
N  ’  N 


2 

N 


N-l 

N 


where  p  is  the  distance  from  the  center  of  the  disk,  a  is  the  radius  of 
the  disk,  and  N  is  a  positive  integer.  BK  is  the  dimensionless  product 
ka  where  k  is  the  wave  number. 

The  only  array  in  the  program  is  CQ.  The  minimum  allocation 
for  CQ  is  given  by 

COMPLEX  CQ(N) 

The  program  prints  out  the  real  part,  the  imaginary  part,  and  the 
magnitude  of  C^  of  (180)  under  the  heading  "REAL  JR  IMAG  JR  MAG  JR." 

The  real  part,  the  imaginary  part,  and  the  magnitude  of  C^  of  (181)  are 
printed  out  under  the  heading  "REAL  JP  IMAG  JP  MAG  JP."  The  real  part, 
the  imaginary  part,  and  the  magnitude  of  C^  of  (190)  are  printed  out 
under  the  heading  "REAL  Q  IMAG  Q  MAG  Q."  The  ith  row  of  numbers  printed 
out  under  any  one  of  these  three  headings  is  for 


p/a  =  (i-l)/N 
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The  sample  input  and  output  data  are  for  the  disk  of  radius 
0.002  wavelengths  with  N  =  30.  The  sample  output  data  for  | C ^ |  are 
plotted  as  the  solid  curve  in  [1,  Fig.  1],  Similarly,  |c^j  is 
plotted  in  [1,  Fig.  2],  and  |C^|  is  plotted  in  [1,  Fig.  3].  The 
solid  curves  were  obtained  by  drawing  straight  line  segments  between 
the  data  at 

p/a  =  0,  1/30,  2/30,  3/30,  ...  29/30 

It  was  possible  to  extend  the  curve  for  |c^|  to  the  rim  of  the  disk 

because  |C  1  is  known  to  be  zero  there. 

1  P1 

The  function  FABS(X)  returns  the  magnitude  of  the  complex 

variable  X.  If  either  of  the  real  and  imaginary  parts  of  X  is  less 
-35 

than  10  in  magnitude,  then  its  contribution  to  the  magnitude  of  X 
is  omitted  in  order  to  avoid  a  machine  underflow. 

Inside  DO  loop  19  in  the  main  program,  Cp  of  (180)  is  accumu¬ 
lated  in  CR,  C,  of  (181)  is  accumulated  in  CP,  and  C  of  (190) 

<p  q 

is  accumulated  in  CQ(J).  The  index  J  of  DO  loop  19  obtains 

p/a  =  (J-D/N 

In  accordance  with  (164),  line  28  puts  p/a  in  R.  Line  31  puts  of 

(183)  in  SR,  line  32  puts  of  (185)  in  SP,  and  line  33  puts  of 

(192)  in  SQ.  Line  34  puts  C  of  (182a)  in  CR,  line  35  puts  C^  of 

-20 

(184a)  in  CP,  and  line  36  puts  C  ^  of  (191a)  in  CQ(J).  If  ka  <  10  , 

nothing  is  added  to  CR,  CP,  and  CQ(J)  for  fear  of  machine  underflows. 
-20 

If  ka  _>  10  ,  lines  40-42  add  C^*  C^?,  anc*  ^q2  to  anc*  ^Q(J) , 

respectively.  If  ka  ^  10  ^ ,  lines  45-47  add  C^,  C^»  anc^  ^q3  to 


CR,  CP,  and  CQ(J).  If  ka  >  LO  ,  lines  50-52  add  C  .,  C, . ,  and  C  , 

p4  <J>4  q4 

to  CR,  CP,  and  C0(J).  If  ka  >  10  lines  55-57  add  C  c,  C,  ,  and  C  c 

—  po  q)5  q5 

to  CR,  CP,  and  CQ(J).  Line  58  puts  the  magnitude  of  CR  in  SR,  and  line 
59  puts  the  magnitude  of  CP  in  SP. 

DO  loop  17  prints  out  CQ(J). 


001 c  LISTING  OF  THE  COMPUTER  PROGRAM  THAT  CALCULATES  B0UWKAMP#S 

002C  POWER  SERIES  SOLUTION  FOR  A  CIRCULAR  DISK 

003 //PGM  JG6  (XXXX.XXXX.1.1 ) .* MAUTZ. JOE*.REGION=200K 

004//  EXEC  WATFIV 

005//GO.SYSIN  OO  4 

006  S  JOB  MAUTZ.T IME=1 .PAGES=40 

OOT  FUNCTION  FABS(X) 

008  COMPLEX  X 

009  R= ABS( REAL (X ) ) 

010  A=ABS( AIMAGIXI I 

Oil  IFIR.LT.  UE-351  R=0. 

012  IFIA.LT.1.E-3S)  A=0. 

013  FABS=SORT IR4RFA* Al 

014  RETURN 

015  END 

016  COMPLEX  U.CR.CP.CCI 180) 

017  C=4./3. 141593 

018  U=<0.*1.) 

019  RE ADI  1 • 1 2 )  N.BK 

020  12  FORMAT! I3.E14.7) 

021  WR I TE (3*14)  N.BK 

022  14  FORMAT  < • 0  N • .6X . • BK • / I  X . 1 3. E 1 4 . 7 ) 

023  0=1 ./N 

024  WRITE (3* 13) 

025  13  FORMAT! *0  REAL  JR  IMAG  JR  MAG  JRa.6X.*REAL  JP  IMAG  JP 

026  1  MAG  JPa) 

027  DO  19  J=1.N 

028  R=<  J— 1 f  40 

029  R2=R4R 

030  S=SQRT< 1.-R2) 

031  SH=C*S 

032  SP=-C/S 

033  SQ=-SP4R 

034  CR= 1 »333333*SR4BK4U 

035  CP=.66666674SP4!2.-R2)48K4U 

036  CQ(J)=S042. 

037  IF<BK  .LT.l .E— 20)  GO  TO  11 

038  BK2=BK*BK 

039  BK3=BK24BK 

040  CR=. ae88889E-01 4SR4< 7. -R2 ) 4BK34U+CR 

041  CP=. 2222222E— 01 *SP4 (P24(R2— 8.1428. )4BK34U4CP 

042  COf J)  =  SQ4BK2M4.-R2)/3.4CQI J) 

043  lFCBK.LT.l.E-15)  GO  TO  11 

044  EK4=BK3*8K 

045  CR=. 2829421 4SR4BK44CR 

046  CP=.28294214SP4BK4+CP 

047  CQIJ)=— •  5658  84  2*  S04BK34U+CC { J9 

048  IF  <  BK .LT.l.E— 12)  GO  TC  11 

049  BK5=BK44BK 

050  CR=-1*/3150.4SP4< ( 188.— 9.4R2) 4R 2—664* ) 4BK54U+CR 

051  CP=-1 ./6300.4SP4 (((3.4R2-54.) 4R2+704. ) 4R2-1328. ) 4BK54U4CP 

|0S2  CQ< J)=SO*BK44(R24l3.4R2-44.)488.)/180.4CQ< J) 

053  IFIBK.LT.  1.E-10I  GO  TC  11 

054  0K6=8K54BK 

055  CR=-. I886231E-024SR4C 15.4R2-148 . )40K64CR 

056  CP=-.  1 88628 IE—  024  SP4f  29.4R2— 148. ) 4BK64CP 

057  CQI J)  =  S04BK 54. 22C3537E-0 l4f-26. +5. 4R2)4U4CQC  J) 

|058  11  SR=FABS  <CR  ) 

059  SPsFABS(CP) 

060  WRI TE (3.15)  CR.SR.CP.SP 


-061  15  FORMAT! 1X.3E1 1.4 

.  062  19  CONTINUE 

. 2X.3E1 1 .4) 

063  MR I TE (3 • 18} 

064  1 8  FORMAT < • 0  REAL 

Q  IMAG 

Q  MAG  Q*> 

*  065  00  17 

J—  1  *N 

r  066  SC— FA8SIC0C  Jl • 

067  MRITE13.15)  CCHJl.SO 

.  068  17  CONTINUE 

069  STOP 

\  070  END 

|  SOATA 

F  30  0 • 1 256637E— 0 1 

SSTOP 

/• 

✓/ 

fl  PRINTEO  OUTPUT 

t*  N  EK 

l  30  0.1 25C637E— 0 1 

REAL  JR 

■MAC  JR 

MAG  JR 

REAL  JP 

IMAG  JP 

MAG  JP 

0.  8985E— 08 

0.2I33E-CI 

0.21 33E-01 

— 0.8985E— 08—0 

•  21 33E-0  1 

0.21 33E— 01 

«  0.8980E-08 

0.2132E-01 

0.2132E-01 

— 0 • 8 99 OE -08—0 

.21 33E-0 1 

0.2133E— 01 

j  0* 8965E—  08 

0.2129E-01 

0.21 29E-01 

-0.9005E-08-0 

•21 33E-0 1 

0.2133E-01 

0.8940E-08 

0.2123E-01 

0.2123E-01 

-0.9030E-08-0 

.2134E-01 

0  .2134E-01 

0.8905E-08 

0.21 14E-01 

0.21 14E-01 

—0. 9066E— 08— 0 

•21 34E— 0 1 

0.2134E-01 

0* 8859E—  C8 

0*21 046—01 

0.21  04E-01 

-0.91 12E-08-0 

•2134E-0I 

0.2134E-01 

0.8803E-08 

0.2090E-01 

0.2090E-01 

-0.91 70E-0e-0 

•2134E-01 

0.2134E— 01 

1  0.8737E-08 

0.2075E-01 

0.2075E-01 

— 0.9 240E— 08-0 

•2134E-01 

0.2134E-01 

*  0. 8660E-08 

0  .2056E-01 

0.2056E-01 

-0.9322E— 08— 0 

. 2 1 35E— 0 1 

0 . 2135E— 0 1 

0.8571E-08 

0.2035E-01 

0.2035E-01 

—0.941 9E-08— 0 

.2136E-01 

0.2136E-01 

0.84  71E-C8 

0*201  IE— 01 

0.20 1 l E— C 1 

—0 .9530E— 08— 0 

•2137E— 01 

0.2137E— 01 

0.8359E-08 

0*  1985E-01 

0. 1985E-01 

— 0.9657E-08— 0 

•2139E-01 

0.2139E-01 

0.  8235E-08 

0. 1955E-01 

G.1955E-01 

— C.9803E— 08— 0 

•2142E-01 

0.214  2E— 0 1 

1  0.8097E-08 

0.1923E-01 

0.I923E-01 

— 0.9970E-08-0 

.2I45E-01 

0.214  5E—  0 1 

*  0*  7947E—08 

0*1 887E—  0 1 

0.188 7E— 01 

-0.1 016E-07-0 

•21 50E-0 1 

0 .2 150E-01 

0.7781E-08 

0.1848E-01 

0.1 84  8E— 0 1 

-0. 1037E-07-0 

•2156E-01 

0.2156E-01 

0.7600E-08 

0.1805E-01 

0  •  1 805E  —  0 1 

— 0.1062E-07— 0 

•21 63E— 0 1 

0.2163E— 01 

0*  7403E-08 

0.1 7S8E-01 

0.1 758E-01 

— 0.1090E-07-0 

.21 74E-01 

0.2174E— 01 

0.7188E-08 

0.1 707E-01 

0.1707E-01 

-0.  1 1 23E— 07— 0 

•2187E-01 

0.2187E-01 

N  0.6953E-08 

0.1 65 IE- 01 

0  •  1 6  E 1 E  — 0  1 

-0. 1 161 E— 0  7—0 

. 2204E-0 1 

0.2204E— 01 

0.6697E-08 

0.I590E-01 

0.1 590E-0 1 

—0.1 205E— 07— 0 

•2226E— 01 

0.2226E-01 

0.6416E-08 

0.1S24E-01 

0.1524E-01 

-0. 1258E-07-0 

• 2256E— 0 1 

0.2256E-01 

0.6109E-08 

0.1 450E-01 

0.1450E-01 

-0. 1322E-07-0 

. 2294E-01 

0.229  4E  —  01 

0.5769E-08 

0.I370E-01 

0. 1370E-01 

-0. 1399E-07-0 

. 2346  E—0 1 

0.2346E-01 

l’  0*  5391 E— 08 

0.1230E-01 

0.12  S0E-01 

-0.1497E-07-0 

.241 8E-01 

0.2418E-01 

'*  0.4967E-08 

0.1179E-01 

0.11 79E-01 

-0. 1625E-07 -0 

.252  OE- 0 1 

0.2520E-01 

0 • 4482E— 08 

0.1064E-01 

0. 1064E-01 

—  0.1801 E-0  7— 0 

•2670E-01 

0 . 2670E— 0 1 

0.3916E-08 

0.9300E— 02 

0.9300E— 02 

—  0.2061 E-0  7-0 

•2912E-01 

0.2912E-01 

0. 3226E—  08 

0.7659E-02 

0. 7659E-02 

— 0.2503E-07— C 

.  3354E-01 

0.3354E-0I 

0.2300E-08 

0.5462E— 02 

0 . 5462E— 0  2 

—0.3 50 9 E-0 7—0 

.  4440E  — 0 1 

0.444  OE— 0 1 

J  REAL  0 

IMAG  Q 

MAG  Q 

O.OOOOE+OO 

O.OOOOE+OO 

O.OOOOE+OO 

0.8494E-01- 

0.4769E-07 

0.8494E-01 

0.1702E+00- 

0.9555E-C7 

0. 1 702E+00 

0.2S60E+00- 

0. I437E-06 

C.2560E+00 

■.  0.3426E+00- 

J  0.4305E +  00- 

0.1924E-06 

0.3426E+00 

0.241 7E-06 

0.43  05E  +  00 

0*51 99E+00- 

0.2919E-06 

0.51 99E+00 

0*611 1 E  ♦00— 0  *343  IE  — 36  C.61UE  +  00 
0*7 046E400  — 0  *3957E  —  06  0.7046E+00 
0*  8009E400— 0  *4497E— 06  0.8009E+00 
0*9004E*00-0*5056E-06  0.9004E+00 
0*1 004E+0I — 0*5636E— 06  0.1004E+01 
0* 1 1 1 1E40I-0.6241E-06  0.1111E+01 
0*  1  225E  +  0I -0.6876E-06  0.1225E+01 
0. 1344E401-C.7545E-06  0.1344E+01 
0. 1470E*01-0 .325eE-06  0.1470E+01 
0. 1606E401-0.9016E-06  C.1606E+01 
0.  1751E+Ol-0*9e35E-06  0.175IE+01 
0. 1910E401-0.I072E-05  0.1910E401 
0* 2084E+0I— 0* 1 1 70E-05  0.2084E+01 
0.2278E401— 0.1275E-05  0.2278E+01 
0.2496E+01-0 .1402E-05  0.2496E+01 
0*2  747E  +  01 -0*1542E  — 05  0.2747E+01 
0.3041E+01-0.1 70eE-05  0.3041E401 
0. 3396E401-0.1907E-05  0.3396E401 
0*  3839E+01 — 0*2156E— 05  0.3S39E+01 
0.4424E+CI-0.2484E-C5  0.4424E+01 
0.5258E401-0.2953E-0S  C.5258E+01 
0*662 IE 4 01 -0*37ieE— 05  0.6621E+01 
0*  9615E+0I— 0*539 9E- 05  0.9615E+01 
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VIII.  THE  MIE  SERIES  SOLUTION  FOR  A  SPHERE 


The  known  solutions  for  the  electric  current  and  electric  charge 
on  a  conducting  sphere  excited  by  an  incident  plane  wave  are  plotted  in 
[1,  Figs.  4,  5,  and  6].  The  known  solution  for  the  electric  current  was 
obtained  from  the  Mie  series  solution  [4,  Eq.  (6-103)].  Actually,  [4, 

Eq.  (6-103)]  had  to  be  modified  so  as  to  be  valid  for  the  incident  plane 
wave  [1,  Eq.  (110)]  which  travels  in  the  negative  z  direction  rather  than 
for  the  incident  plane  wave  [4,  Eq.  (6-96)]  which  travels  in  the  positive 
z  direction.  The  electric  charge  was  obtained  via  the  equation  of  con¬ 
tinuity.  The  foregoing  electric  current  and  electric  charge  were  calcu¬ 
lated  by  means  of  a  computer  program.  This  program  is  listed  at  the  end 
of  this  section. 

If  Eq  is  set  equal  to  n>  if  P^(cos  0)  is  replaced  by  -P^(cos  0), 

1 '  l' 

and  if  P  (cos  0)  is  replaced  by  -P  (cos  0),  then  [4,  Eq.  (6-103)] 
n  n 

becomes 


where 


J(0,<{>)  =  UqJq(0)  cos  4>  +  u^J^(0)  sin  <p 


t  ra\  1  r  d8  n 

n=l 


'0""  ka  J;,  “n  '  fi(2)' 


j  ~  P*(cos  0)  P^cos  0) 
”*  _  ^  n 


(ka) 


V°>  ■  S  an  (- 

r  n=l 


-jP*(cos  0) 


sin  0  H(2)(ka) 
n 

35  pi<cos  0> 


) 


sin  0  fi(2)’(ka)  fi(2)(ka) 

n  n 


) 


(193) 


(194) 


(195) 


where 


a  =  j ~n (2n+l) 
n  n(n+l) 


(196) 


In  (193),  u^  and  are  the  unit  vectors  in  the  0  and  <p  directions,  respec¬ 
tively.  As  given  by  (193),  .1(0, 4>)  is  the  electric  current  on  the  conducting 


sphere  of  radius  a  centered  at  the  origin  and  excited  by  the  plane 

~inc 

whose  electric  field  is  _E  given  by 


pine  -jkz 

E  =  u  n  e 

—  — x 


(197) 


1  1  1 '  1 ' 

It  was  necessary  to  replace  Pn  by  -P^  and  P^  by  -P  in  [4,  Eq.  (6-103)] 

because  the  P^  used  in  the  present  report  is  taken  to  be  the  negative  of 
n 

the  P^  used  in  [4],  The  P1  used  in  the  present  report  is  that  which  is 
n  n 

calculated  in  [5].  The  right-hand  sides  of  (194)  and  (195)  coincide  with  the 
right-hand  sides  of  [5,  Eqs.  (2)  and  (3)]. 

The  electric  charge  associated  with  the  electric  current  (193)  is 
called  qe(0,<J>).  As  calculated  from  [1,  Eq.  (1)]  with  V  •  given  by 


[1,  Eq.  (B-3) ] , 


q  (0  »d>)  =  cos 


(198) 


where 


^(9)  =  k^Tiire  (ie  (sin  9  V9))  '*■  Ve» 


Substitution  of  (194)  and  (195)  into  (199)  produces 


(199) 


i (0) - hr  l 


oo  a 


P  (cos  0) 


rrz  \  S2T^(-^^<sln9^p^cose>)+Tir">  <200) 

(ka)  n=l  H  (ka)  sin  0 

n 


Thanks  to  the  associated  Legendre  equation  [4,  Eq.  (E-l)] 


■  v1— a  (sin0-^  P1(cos  0))  +  [n(n+l)  -  — 1  P^(cos  0)  =  0 

sin  0  d0  d0  n  .  2  n 

sin  u 


(201) 


expression  (200)  reduces  to 


,  «=  a  n(n+l)P  (cos  0) 

^  \  I  v  n  _ n _ 

q(0)  -  2  L  (2)’ 

(ka)  n=l  ft;  (ka) 


(202) 
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So  far,  it  has  been  established  that  the  incident  plane  wave  (197), 
which  travels  in  the  positive  z  direction,  induces  the  electric  current  (193) 
and  the  electric  charge  (198)  on  the  sphere.  However,  the  objective  is  to 
find  the  electric  current  and  the  electric  charge  induced  on  the  sphere  by 
an  incident  plane  wave  that  travels  in  the  negative  z  direction. 

If  the  situation  in  which  the  electric  current  (193)  exists  on  the 
sphere  excited  by  the  incident  electric  field  (197)  is  rotated  180°  about 
the  x  axis,  then  the  resulting  situation  is  that  of  the  electric  current 

induced  on  the  sphere  excited  by  the  incident  plane  wave  whose  electric  field 

.  _inc  .  , 

is  E  given  by 


inc  ikz 

E  =  u  H  eJ 
—  —x 


(203) 


This  plane  wave  travels  in  the  negative  z  direction  and  coincides  with 
(1,  Eq.  (110)1.  Therefore,  the  electric  current  J  induced  on  the  sphere 
in  [1]  is  £  of  (193)  rotated  180°  about  the  x  axis. 


Substitution  of 


u«  =  u  cos  0  cos  <b  +  u  cos  9  sin  <)>  -  u  sin  9 
~ x  ~ y  — z 


u,  =  -  u  sin  <p  +  u  cos  d> 
— $  — x  — y 


into  (193)  gives 


(204) 

(205) 


J(0,(|>)  =  (u^cos  9  cos  <j>  +  u^cos  9  sin  <j>  -  u^sin  0)  Jq(9)  cos  <J> 


+  (-u^sin  4>  +  videos  <{>)  J^(0)  <P  (206) 

Expression  (206)  is  more  suitable  than  (193)  for  180°  rotation  about  the 
x  axis.  A  180°  rotation  about  the  x  axis  amounts  to  replacement  of  x 
by  itself,  replacement  of  y  by  -y,  and  replacement  of  z  by  -z.  This  is 


written  as 
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x  -*■  x 


y  -y 


z  -*■  -z 


(207) 


Now,  (207)  is  equivalent  to 


0  ■>  TT  -  0^1 

4>  -  <j>  j 


(208) 


When  rotated  180°  about  the  x  axis,  u  remains  unchanged,  u  changes  to 

x  -y 

A 

-u  ,  and  u  changes  to  -u  .  It  is  now  evident  that  J[(0,4>)  of  (206)  can 
y  z  z 

be  rotated  180°  about  the  x  axis  by  taking  the  value  of  J(0,<t>)  with  the 
signs  of  its  y  and  z  components  changed  and  placing  this  value  not  at 
(0,<J>)  but  at  (tt-0,  — <J>).  Hence,  if  the  result  of  180°  rotation  of  £(0,<j>) 
about  the  x  axis  is  called  J(0,<f>),  then 


J_(tt-0,  -4>)  =  (^cos  0  cos  4>  -u^cos  0  sin  <J>  +  u^sin  0)  Jg  (0)  cos  cj> 


+  (-u  sin  (b  -  u  cos  (b)  J.(0)  sin  d> 
— x  — y  <p 


(209) 


Replacement  of  0  by  tt-0  and  4>  by  -<p  in  the  arguments  of  J^,  Jg,  J^,  the 
sines,  and  the  cosines  in  (209)  yields 


J(0,<J>)  =  -  (u  cos  0  cos  (p  +  u  cos0sin4>-u  sin  0)  Jfl(iT-0)  cos  4> 


/A 

+  (-u^sin<}>+  HyCOS  $)  s^-n  ^ 


(210) 


Now,  (204)  and  (205)  reduce  (210)  to 


,J(0,(j>)  =  -  UgJg(TT-8)  cos  $  ■*"  (tt— 9 )  sin  <}> 


(211) 


The  electric  charge  associated  with  the  electric  current  (211)  is 


called  q  (0  ,<}>) .  As  calculated  from  [1,  Eq.  (1)]  with  V  •  given  by 

6  S 

H,  Eq.  (B-3)  ] , 
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I 


> 


i 


4 


*  . 


qe(6 ,4>)  =  q-(-TTc9)  cos  <p  (212) 

where  q(0)  is  given  by  (199)  which  was  previously  reduced  to  (202). 

According  to  the  development  (20J)-(212),  the  incident  plane  wave 
(203),  which  travels  in  the  negative  z  direction,  induces  the  electric 
current  (211)  and  the  electric  charge  (212)  on  the  sphere.  The  computer 
program  listed  at  the  end  of  this  section  calculates  the  functions  Jg(0), 
J^(0),  and  q (0)  that  appear  in  (211)  and  (212).  Of  course,  these  functions 
also  appear  in  expressions  (193)  and  (198)  for  the  electric  current  and 
electric  charge  induced  on  the  sphere  by  the  incident  plane  wave  (197),  which 
travels  in  the  positive  z  direction. 

The  computer  program  listed  at  the  end  of  this  section  contains  the 
subroutines  BES,  LEG  and  CURNTC.  The  subroutines  BES  and  LEG  are  exactly 
the  same  as  in  [5]  and  are  described  in  detail  there. 

The  subroutine  CURNTC (X,  NT,  CUR,  Q)  puts  JQ(0)  of  (194)  for  0  = 
(K-1)tt/(NT-1)  in  CUR(K)  for  K-l, 2, ..  .NT-  Similarly,  CURNTC  puts  5^(0)  of  (195) 
in  CUR(NT+1)  through  CUR(2*NT)  and  q(0)  of  (202)  in  Q(l)  through  Q(NT) . 

The  value  of  ka  in  (194),  (195),  and  (202)  is  equal  to  X.  The  arguments 
X  and  NT  are  input  arguments  and  the  arguments  CUR  and  Q  are  output  argu¬ 
ments.  The  subroutine  CURNTC  calls  the  subroutines  BES  and  LEG.  Minimum 
allocations  in  the  subroutine  CURNTC  are  given  by 

COMPLEX  CUR(2*NT) ,  Q (NT ) ,  H(N) ,  HP(N),  HQ(N) 

DIMENSION  BJ(2*N+1),  BJP(2*N+1),  BY(N+1),  BYP(N+1) 

where 

N  =  MIN(5.+4.*X,  29.)  (213) 

Here,  MIN  denotes  minimum  value. 
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The  subroutine  CURNTC(X,  NT,  CUR,  Q)  was  created  by  appending  Q 
to  the  list  of  arguments  of  the  subroutine  CURNTC(X,  NT,  CUR)  listed  in 
[5,  p.  11]  and  calculating  Q  by  means  of  (202)  .  For  the  calculations,  the 
oo  N 

l  in  (194),  (195),  and  (202)  was  replaced  by  £  where  N  is  given  by  (213). 
n=l  n=l 

In  the  subroutine  CURNTC,  statement  14  puts  the  spherical  Bessel 

function  of  the  first  kind  j  (ka)  in  BJ(n+l),  j ' (ka)  in  BJP(n+l),  the 

n  n 

spherical  Bessel  function  of  the  second  kind  y^(ka)  in  BY(n+l),  and 
y^(ka)  in  BYP(n+l)  for  n  =  0,  1,2,...N.  The  functions  ’  (ka)  and  (ka) 

in  (194),  (195),  and  (202)  are  given  in  terms  of  the  spherical  Bessel  func¬ 
tions  and  their  derivatives  by 

H^(ka)  =  ka(j  (ka)  -  jy(ka))  (214) 

n  n 

(ka)  *  kaj ' (ka)  +  j  (ka)  -  j  (ka  y'(ka)  +  y  (ka))  (215) 

n  n  n  n  n 

a 

DO  loop  11  puts  - -  in  HC-.)  for  n  =  J.  Furthermore,  DO  loop  11 

ka  H'2 ' (ka) 
n 

ja  a  n(n+l) 

puts  - -  in  HP(J)  and  - 5 — ppr~, -  in  HQ(J).  The  index  K 

ka  H(2)  (ka)  (ka)2  H(2)  (ka) 

n  n 

of  DO  loop  12  obtains  9  according  to  0  =  (K-l)71/ (NT-1).  Statement  15  puts 

P^(cos  9)  in  BJ(N+l+n)  and  P^(cos  0)  in  BJP(N+l+n)  for  n  *  1,2,...N. 
n  do  n 

DO  loop  13  accumulates  Jq(0)  of  (194),  .5^(0)  of  (195),  and  q(0)  of  (202) 
in  Gl,  G2,  and  G3,  respectively. 

The  main  program  reads  input  data  from  a  punched  card  according 
to 

READ(1,23)  NT,  XC 
23  FORMAT ( 13 ,  F14.7) 
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For  ka  =  XC,  the  main  program  prints  out  the  quantities  Jg  (0)  of  (194), 

J.(0)  of  (195),  and  q(0)  of  (202)  at 

0  =  (i-l)Tr/  (NT-1)  ,  i  =  1,2,... NT  (216) 

where  NT  is  a  positive  integer  and  NT  >  1.  The  main  program  requires  the 
subroutines  CURNTC,  BES,  and  LEG.  Minimum  allocations  in  the  main  program 
are  given  by 

COMPLEX  CUR(2*NT) ,  Q(NT) 

The  main  program  prints  out  the  real  part  and  the  imaginary  part 
of  3Q(0)  under  the  heading  "REAL  JT  IMAG  JT."  The  real  and  imaginary 
parts  of  3^(0)  are  printed  out  under  the  heading  "REAL  JP  IMAG  JP."  The 
magnitude  of  JQ(0)  is  printed  out  under  the  heading  "MAG  JT."  The  magni- 
tude  of  J,(0)  is  printed  out  under  the  heading  "MAG  JP."  Down  farther, 

<f> 

the  real  part,  the  imaginary  part,  and  the  magnitude  of  q(0)  are  printed 
out  under  the  heading  "REAL  Q  IMAG  Q  MAG  Q."  The  ith  row  of  numbers 
printed  out  under  any  one  of  these  headings  is  for  0  given  by  (216).  The 
sample  input  and  output  data  are  for  the  sphere  of  radius  0.002  wavelengths 
with  NT  =  31.  The  sample  output  data  for  |.Iq(0)|  are  plotted  as  the  solid 
curve  in  [1,  Fig.  4].  The  solid  curve  was  obtained  by  drawing  straight  line 
segments  between  the  data  at  the  values  of  0  given  by  (216).  Actually, 

|  Jg  (0)|  was  plotted  at  t  =  0  in  [1,  Fig.  4],  and  I  Jg  (tt)  |  was  plotted  at 
t  =  Tra  in  [1,  Fig.  4].  However  t  =  0  corresponds  to  0  =  it,  and  t  =  ira 
corresponds  to  0  =  0  so  that  the  solid  curve  in  [1,  Fig.  4]  is  really  a 
plot  of  )  Jg  (tt— 0)  |  .  This  is  in  harmony  with  expression  (211)  for  the 
electric  current  on  the  sphere  excited,  as  in  [1,  Fig.  4],  by  the  incident 
plane  wave  (203),  which  travels  in  the  negative  z  direction.  Similarly, 
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the  sample  output  data  for  | J  (tt— 0 )  |  are  plotted  as  the  solid  curve  in 
[1,  Fig.  5],  and  the  sample  output  data  for  | q (7T— 0 )  |  are  plotted  as  the 
solid  curve  in  [1,  Fig.  6]. 

If  ka  £  0.001,  then  the  subroutine  CURNTC  can  not  be  used  to 

calculate  Jq(6),  J^(9),  and  q(0)  because  the  subroutine  BES  called  by  it 

does  not  calculate  any  of  the  functions  j ' (ka) ,  y  (ka)  and  y'(ka)  whenever 

n  n  n 

ka  £  0.001.  If  ka  is  small,  then  the  spherical  Bessel  functions  j^Cka) 
and  y^(ka)  can  be  approximated  by  (6,  Eqs.  (10.1.4)  and  (10.1.5)] 


j  (ka)  =  - —  (kg) 
n  1*35...  (2n+l) 


x  -1-3-5.  ..(2n-l) 

y  (ka)  - - 7^ - — 

n  ,.  Nn+1 

(ka) 


(217) 

(218) 


Therefore,  when  ka  is  small,  H(2)(ka)  of  (214)  and  H(2)’(ka)  of  (215) 

n  n 


can  be  approximated  by 


g(2)(ka)  =  i  (1-3-5. .. (2n-l) ) 
n  (ka)n 


fi(2)  '  (ka)  =  rjnU-3-5  (2n-l)l 
n  (ka)n+1 


(219) 

(220) 


When  ka  is 
in  each  of 
and  (202). 
defined  by 


In  view  of 


small,  it  suffices  to  retain  only  the  term  for  which  n  ■  1 
the  infinite  series  on  the  right-hand  sides  of  (194),  (195), 
The  required  P^(cos  9),  being  the  negative  of  the  P^(cos  0) 
[4,  Eq.  (E-17)],  is  given  by 

P^(cos  0)  =  t[\  -  cos20  =  sin  (221) 

(219),  (220),  and  (221),  retention  of  only  the  term  for  which 


n  =  1  in  each  of  the  infinite  series  on  the  right-hand  sides  of  (194), 
(195),  and  (202)  gives 
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JQ(0)  =  -1.5 

(222) 

J.(0)  =  1.5  cos  0 

<P 

(223) 

q (0)  =  3  sin  0 

(224) 

If  ka  £  0.001,  then  the  branch  statement  on  line  198  in  the  main 
program  causes  DO  loop  16  to  be  executed.  The  index  J  of  DO  loop  16 

obtains  0  =  (J-I)tt/ (NT-1) .  Inside  DO  loop  16,  JQ(0)  of  (222)  is  put  in 

o 

CUR(J),  J,(0)  of  (223)  is  put  in  CUR(.I+NT),  and  q(0)  of  (224)  is  put  in 
<P 

Q(J). 


If  ka  >  0.001,  then  the  branch  statement  on  line  198  causes 
statement  15  to  be  executed.  For  the  NT  values  of  0  given  by  (216), 
statement  15  puts  J0(0)  of, (194)  in  CUR(l)  through  CUR(NT),  J^(0)  of 
(195)  in  CUR(NT+1)  through  CUR(2*NT),  and  q(0)  of  (202)  in  Q(l)  through 
Q(NT).  DO  loop  10  prints  out  J0(0)  and  .1^(0) .  DO  loop  13  prints 
out  3(0). 


001  C  LISTING  OF  THE  COMPUTER  PROGRAM  FOR  THE  MIE  SERIES  SOLUTION 

:  002  c  FOR  A  SPHERE 

V'  003 //PGM  JOB  ( XXXX.XXXX. 1 • 1 ) .• MAUTZ. JCE« ,REGION=200K 
004  //  EXEC  4ATFIV 
If  005  //GO-SYSIN  DO  * 

006  * JOB  MAUTZ.TIME=1 »PAGES=40 

007  C 

008  C  LISTING  OF  THE  SUBROUTINE  BES 

009  SUBROUTINE  BES IL  *LD. B JP  »  BY« BYP ) 

010  DIMENSION  BJiSO) . 6 JP( 50 ) . B Y( 50 ) . B YP( SO) • BS( 50 > 

■  Oil  L1=(L-1)*NJ 

B012  L3=(LD-ll*NJ 

013  6  IF(XJ-l.E-3)  3.3*4 

014  3  Jl=Ll4l 


015 

J2=L14NJ 

016 

DO  5  J=J1.J2 

017 

8J< J)=0. 

018 

s 

CONTINUE 

019 

BJ« J1 )=1. 

020 

RETURN 

021 

4 

SN=SI  N<  XJ) 

022 

CS=COS(  XJ1 

023 

IFIXJ-15.I  11.12.12 

024 

12 

B J IL 1 ♦ 1 >=SN/XJ 

025 

B J (LI 42 )=(BJ(L141 l-CSI/XJ 

026 

DO  14  1=3. NJ 

027 

I 3=L 1 4 I 

028 

12=13-1 

029 

11=13-2 

030 

B J I  13 1  =  FLO AT (2*1 —3 ) 7X J4B Jf 121-BJf 

031 

14 

CONTINUE 

032 

B3=FL0AT( 2*NJ-1 1/X J*8J( I3I-BJI 12) 

033 

GO  TO  15 

034 

11 

NBJ=XJ*22« 

035 

IF( NBJ.LE.NJ )  NB  J=N  J4 1 

036 

BS( NB J )=0. 

037 

BSINBJ-l )=l.E— 50 

038 

NB J2=NB J— 2 

039 

00  193  1=1. NB J2 

040 

I  2=  N  B  J—  I 

041 

13=1241 

042 

I  1=12-1 

043 

FI=FL0AT(2*I14l)/XJ 

044 

BS( I 1 »=BS( I2)*FI— BSI 13) 

045 

193 

CONTINUE 

046 

Bl=SN/XJ 

047 

B2=E1 /X J-CS/X J 

048 

IF(A8S(Bl )-ABS(B2)> 1.2.2 

049 

2 

BB=B1/BS( 1) 

050 

GO  TO  9 

051 

1 

BB=B2/BS( 2 ) 

052 

9 

00  194  1=1. NJ 

053 

I 1=L1 4! 

054 

B J ( 1 1 )=BS( I )*BB 

055 

1 54 

CONTINUE 

056 

B3=BS( N  J  41) 4B8 

057 

15 

BY(L1 41 )=-CS/XJ 

058 

BY(L142)=(BY(L141 )-SN)/XJ 

059 

00  64  I=3.NJ 

AAA 

I3=L14I 
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2=13-1 

1=13-2 

r< 13  l=FLOATI  2*1 -3>/XJ*BV( I2)-BY( 111 
3NTINUE 

1=FL0  AT  C  2*N J  — 1 >/X J*BV< I3)-BY< 121 

-(ID.EQ.2)  RETURN 

J1  =N  J—  I 

l=L3*l 

2=1.1  *2 

>P<JI)=— BJ< J2> 

rPI J1 )=-BYI  J2 1 

3  65  J=2.NJ1 

>=L3*J 

l=Ll ♦ J— 1 

)=J1*2 

J=2*  (  2  *  J— 1  I 

JP1 J2)=.5*IBJIJ I )  — BJ  I  J3  1 )-CBJ< J1)*BJ< J3IJ/FJ 
tPi J2)=.5*(BY1 J1 I —BY  < J3))-(BY<J1)*BYCJ3)1/FJ 
3NTINUE 
J=F  J*4  • 

2=J2*1 
l  —  J 1  *1 

IP! J2)=.5*IBJ{ J1 ) — B3 )— IBJIJ1 )+B3)/FJ 
rP  (  J2)  =  .  5*  (BY!  J  I  I-B4  )-(BV( J1 >*B4)/FJ 
'TURN 
40 

[STING  OF  THE  SUBROUTINE  LEG 
JBKOUTIfsE  LEGILtLOi  ID.NJ.F.XP.P.  PPI 
t  MENS I CN  PCI  81* PI  59) .PP 1 59  I . PS  1 50 > 

:<  i)=i. 

l  =  M*l 

)  7  J=1.M1 

I! J*1 )=PCI J)*FLCAT( 2*J-I I 
JNTINUE 

>=M*NJ-M*<M-Il/2 
f=IL-l  1*L5 
i=(LC-l  )*L5 
>=ABS<  l.-XP*XP) 

=S0RT|X2) 

)  3  J=  1  •  M 1 

!=L2  +  1  J-l ) *NJ-<  J-2)*I J-t 1/2 
1=1  . 

IJ.NE.l)  X3  =  X 1  * • I J— 1 1 
;<  i»=pci  J)*x3 
i I  2 )  =  PC I J* I ) *  XP4X3 
(J.EQ.MI)  GO  TC  14 
M2*  1I=PS( 1 ) 

M2*2)=PS<2 J 
H=N  J- J*  1 
i  4  1=3. NJ1 

*1-2 
=  1-1 

111  =  2. *XP*PSCI2)-PS(I1)  *FLO  AT  (  2  *  J-3  ) /FLO  AT  C  12)  MXP*PS{  I2  3-PSI 11 1 

I J.E0.M1 I  GO  TC  4 
*M2*  I 
J  2  I *  PS I  1 1 
* f I NUE 
I MOE 


) 


4 


> 


I 


1 


4 
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121 
122 
.  123 

i24 

IS  125 

V  126 
127 
:  126 

129 

130 
m  131 
™  132 

133 

134 

135 

136 
'm  137 
™  138 

139 
$  140 

141 

142 

hi  143 

144 

145 

146 
1 47C 
148C 

S,“,c 

150 

151 

152 

153 

154 
■  «S5 

156 
V-  157 
-  158 
I  159 
160 
"4  161 
.  162 

163 

164 
;  165 

166 
"4  167 
168 

169 

170 

171 

172 
*173 

174 

175 

176 

177 
!78 

*179 

-180 

r . 


IF! ID  *EQ*  2)  RETURN 
00  5  J=l.M 

M2=L4 ♦! J-l 1*NJ-I J-2)*l J-l )/2 

M3=M2+L2-L4 

NJI=NJ- J41 

OC  6  I-2.NJ1 

32=1*241 

JI=M34l-NJ4J-l 

J3=M34l4NJ-J 

1FI J.NE.l . AND. J.KE.M)  GO  TO  8 
IFIJ.NE.1)  GO  TO  12 
PPI J2)=-PI J3) 

GO  TO  6 

12  PP<  J2)=.5*IFLOATI  I*t2*J*I-3)  >*P1  J1  1-PS<  1-1)1 
GO  TO  6 

8  PP( J2 )= .5*1 FLOAT ( 1*1 2*J4 1—31 ) *Pt  J 1)-Pf J3) ) 

6  CONTINUE 
J2=M2*1 
J1=M3-NJ+J 
IF(J.NE.i)  GO  TO  13 
PP< J2I=0. 

GO  TO  5 

13  PPI J2)=.5*FLOAT< 2*J-2)*P(J1J 
5  CONTINUE 

RETURN 

END 

LISTING  OF  THE  SUBROUTINE  CURNTC 

THE  SUBROUTINE  CURNTC  CALLS  V- £  SUBROUTINES  BES  AND  LEG 
SUBROUTINE  CURNTC I  X . NT .CUR . 0 ) 

COMPLEX  CURI 146) . Q ( 73 ) • U. G 1. G2. G3. G4. Hf 29 ) . HP( 29) • HOI  29) 
DIMENSION  H J ( 59)  .8 JP< 59 ) . BY I  SO ) • BYP ( 50 ) 
N=MINII5.+4.4X.29.) 

NP=N41 

14  CALL  BES I  1 • 1* 1 »NP« X • B  J . 0 JP.BY. BYP ) 

U=(0..1 .) 

G 1=  1  a 

OO  11  J=1'N 
G1=-G1*U 

G3=FLQAT 12* J+ 1 )/X*Gt 
Jl= J*l 

G2=l. /FLOAT! J*JI ) 463 

G4=l./( X*B JPI J1 ) *B  Jl J 1  )— U  4 1  X  * BYP I J 1 ) ABY I J 1 ) ) ) 

HI J )=G2/f X*(9J( J| )-U*BTI Jl )) ) 

HPI J)=G4#G2*U 
HCl J)=G4/X*G3 
11  CONTINUE 

DT=3. 1  A  1593/1 NT- 1) 

OC  12  K=1 • NT 
Z=COS<  <K-I)*0T) 

16  IFIABSIZ) • GE • 1 •  )  Z= S 1 GNl . 99998. Z ) 

SN=1./SGRTI 1 •  — Z*  Z ) 

15  CALL  LEGI 1« 1 *1 «NP«2.Z «B J.BJP) 

Gi*0. 

G2=0. 

G3-0. 

00  13  J=  1  »  N 
J1»NP* J 
SB=SN*0J( Jl  ) 

G  1=G  1  4BJPI Jl )*HP( J)*SB*HI  J) 


n 


S 


r 


t> 


■  1 

■  A 
1 
1 

1 


I 


I 


» 


i: 
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tat 

182 

183 

184 
!  185 

I  186 

;  187 
-  188 
;  i89 
-  190C 
I  191 C 
192C 

193 

194 

195 

196 
I  197 

198 

199 

200 
201 
202 

)  203 

204 

205 

206 

207 

208 
1 209 

210 
'21 1 
212 
213 
.214 
*215 
216 

217 

218 
219 
-220 


G2«G2-SB4HP< J)-BJP< J1 )4M< Jl 
G3=G3+HC< J)*EJ< J1 ) 

13  CONTINUE 
CUR1K  )=Gl 
CUR!K»NT)=G2 
Q!K)=G3 

12  CONTINUE 
RETURN 
END 

LISTING  OF  THE  MAIN  PROGRAM 
THE  SUBROUTINES  CURNTC.  BES.  ANO  LEG  ARE  REQUIRED 
COMPLEX  CURI292) *Q< 1461.C2 
READ! 1 .  23 )  NT.XC 

23  FORMAT! I3.E14. 7) 

•  R I TE  <3*24 )  NT.XC 

24  FORMAT!*  NT • ,6X • • XC • ✓ IX • I 3.E l 4.7 ) 

IFIXC-1.E-3)  14*14.15 

14  DT=3.141593/!NT-1) 

00  16  J= 1 • NT 
TH= ! J-t »*0T 
CUR! J)=-l .5 
CUP!J+NT)=1.S*C0S!TH» 

Q! J)=3.*SIN!TH) 

16  CONTINUE 
GO  TO  17 

15  CALL  CURNTC !XC*NT  *CUR*Q I 

17  MR  I TE !3  *27 ) 

27  FORMAT! *0  REAL  JT  I MAG  JT  REAL  JP 

1MAG  JP*I 
DO  10  J=1,NT 
P 1=CABS  <  CUR  <  J ) ) 

C2=CUR< J+NT1 
P2=CABS!C2I 

WRITE13.il)  CUR! J) .C2.P1.P2 

11  FORMAT < 1X.6E1 I .4 ) 


1M  AG 


MAG  JT 


i  11  FORMAT < 1X*6E1 I .4 ) 

10  CONTINUE 
I  4RITE13.12) 

>  12  FORMAT! *0  REAL  Q 

DO  13  J=1*NT 
P1=CABS!Q! Jll 
WRITE<3*11)  Cl J ) *P1 
13  CONTINUE 
STOP 
ENO 

•OATA 

31  0.1 25663 7E-0 1 
•  STOP 

/• 

✓/ 


IMAG 


PRINTEO  OUTPUT 
NT  XC 

31  0.1 25663 7E-01 


REAL  JT 
-0.  1500E401 
-0. 1500E401 
-0. 1500E+01 


IMAG  JT 
0.2932E-01 
0 . 29 1 6E— 0 1 
0 .2868E- 0 1 


REAL  JP  IMAG  JP  MAG  JT  MAG  JP 
0. 1500E401-0.2932E-01  0.1500E401  O.I500E401 
0.1492E+01— 0. 2909E—0 1  0.1500E+01  0.1492E+01 
0.1467E+01— 0.2842E-01  0.1500E+01  0.1467E+01 


-0.1500E+01  0.2789E-01  0. 1 426E+01-0.2732E-0 1  0.1500E+01  0.1427Ef01 
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-0.  1500E+01  0.2679E-01  0.  1  3 70E+  0  1-0 . 2586E-0 X  0. 1500E+01  0.1370E+01 
-0.15C0E+01  0.2539E-01  0 . 1 299E >0 1 -0 .240 9E-0 I  0.1500E+01  0.1299E+01 

-0.1500E+01  0  .  23  72E—  0  1  0*  1  2  i  3E  +  0 1-0  •  2209E-9 1  0.1500E+01  0.1214E+01 
-0.  1500E+01  0.2I79E-01  0.  1 1  1 5E  +  0  1  -0.  1  994E-0 1  0.  1500E  +  01  0.11156+01 

—  0*1 500E  + 01  0.1962E-0I  0 . 1 0 04E  +  0 1-0 . 1  7 76E-0 1  0. ISOOE+OI  0.1004E+01 

—  0*  1500E+01  0.I723E-01  0. 88  16E  +  00-0 • 1 56 1E-0 1  0.1500E+01  0.8817E+00 

-0.  1500E+01  0.1466E-01  0.7499E +  00-0 .1361E-0I  0. ISOOE+OI  0.7501E+00 

—  0* ISOOE+OI  0*11 93E- 0 1  0. 6 i 0 1 £400-0 • 1 1 84E-01  0.1500E+01  0.6102E+00 
-0.  1500E+01  0.9060E-02  0 . 46 35E+00-0* i 9 3eE-0t  C.1500E+01  0.4636E+00 
-0.  1500E  +  01  0.6096E-02  C .3 1 1 8E +0 0-0 .9284E-02  O.ISOOE+Ol  0.3120E+00 
—0*1 500E+0 I  0.3064E-02  0. 1 568E+00-0 . 8608E-02  0.1500E+01  0.1570E+00 
—0*1 500E+C1— 0.9828E-06  0 . 4459E-06-0 • 8379E-02  0.1500E+01  0.8379E-02 
-0* ISOOE+OI — 0 • 3066E-02  — 0  *1 568E+00— 0.  8608E  — 02  0.1500E+0  1  0. IS70E+00 
-0* ISOOE+OI  — 0* 609 8E— 02— C*  31  I 8E+00-0 .9285E-02  0.1500E+01  0.3120E+00 
-0* ISOOE  +  OI  — 0  *9062E-02-0 • 4  6  35E  +  00— 0 • 1 038E— 0 1  0. ISOOE+OI  0.4636E+00 
-0.  1500E  +  01  -0*  H93E-01  -C.61  ClE  +  00-0.1 184E-01  0* ISOOE  +  OI  0.6I02E+00 
-0* ISOOE  +  OI  — 0.1 466E- 01-0. 7499E  +  00—0 .1362E-0!  0.I500E  +  01  0.7501E+00 
-0. 15 OOE+O 1-0. 1724E-01- 0.881 6E+ 00-0. 1562E-0 l  0. ISOOE+OI  0.881 7E+00 
-0. 1 500E  +  01  — 0  *  1 962E-0 1-0.1004E  +  0  1-0.1776E-01  0. ISOOE  +  OI  0.1004E+01 
-0. 1S00E  +  01-0.21 79E-0I-0. 1 1 1SE  +  0  1-0. 199SE-01  0. 1500E  +  01  0.1 1 15E  +  01 
-0. 1500E+ 01-0. 2372E- 01-0 . 1 2 1 3E+ 0 1-0.22 09E-0 1  0.1500E+01  0.1214E+01 
-0. ISOOE+OI -0.254 OE- 01-0. 1 299E+0 1-0.2409E-01  0. 1500E+0 1  0.1299E+01 
-0. 1500E+01— 0.26 79E-0 1-0.1 370E+0 1-0.2586E-01  0. ISOOE+OI  0. 1370E+01 
-0. ISOOE+OI -0 .2 789E- 01 -0 . 1 426E+0 1 -0 .2 732E-0 I  0.1500E+0 1  0. 1427E+01 
-0.1500E+Cl-0.286eE- 01-0. 14676+0 1-0.2842E-01  0. ISOOE+OI  0.1467E+01 
-0. 1 500E  +  01  — 0.2916E— 01— 0. 1492E+0  1-0.2909E-0I  C.1500E  +  01  0. 1492E+0I 

—  0.1500E  +  01  — 0. 29 3 2 E—  01 —0 • 1 5 OOE+O 1—0 .2932E— 0 1  0.1500E+0 1  0. 1500E+01 


REAL  Q 
0.  1  899 E— 01 
0.3136E+00 
0.6238E+00 
0. 92716+00 
0. 1220E+01 
0. 1 500E+01 
0. 1 763E+01 
0. 2008E+01 
0. 2230E+01 
0.2427E+01 
0.2598E+01 
0.2741E+01 
0. 2853E+01 
0.2935E+01 
0 . 2984E+  0 1  ■ 
0. 3C00E+01 
0.2984E+01 
0.2935E+01 
0.2653E+01 
0.27416+01 
0. 2558E+01 
0. 2427F+01 
0. 2230E+01 
0.2008E+01 
0. 1 763E+01 
0. 1500E+01 
0. 1220E+01 
0.9271  E  +  00 
0.6238E+00 
0.31366*00 
0. 1 C36E- 02 


IMAC  O 
— 0 . 1 988E— 03 
1—0  •  3266E-  02 
>—  3.6390E  —  02 
>  — 0  •  92  3  4  E— 02 
-0.1 167E-01 
-0.1  361E-01 
-0. 14946-01 
—0 • 1562E-0 1 
—0.1 562E— 0 1 
-0 .1494E-01 
-0.13616-01 
-C.  11686-01 
-0.92376-02 
-0 .6393E-0? 
-0. 32 70E-02 
-0 .3979E-C5 
0 . 326  2  E -02 
0 .6385E-02 
0.9229E— 02 
0.1 167E-01 
0 .1360E-01 
0.1494E-01 
0- 1562E-01 
0. 1562E-0t 
0. 149 4 E  —  01 
0 • ! 360E-01 
0.1 167E-01 
0 .92326-02 
0.6388E-02 
0.3265E-02 
0.1 085E-C4 


MAG  a 
0.1 899E— 0 1 
0.3136E+00 
0. 62386+00 
0 .92  7 1 £+00 
0. 1220E+01 
0.1 5C0E+01 
0.  1764E  +  01 
0.2008E+01 
0 .22  306+0 1 
0 .2427E  +  01 
0.25986+01 
0.27416+01 
0. 28S3E+0 1 
0 • 2935E  +0 1 
0 • 2984E  +0 1 
0. 3000E+0 1 
0.29846+01 
0 .29356+0 1 
0.2853E+01 
0.2741E+01 
0  .2598E+01 
0.2427E+01 
0 • 22J0E+0 1 
0 • 2008E  +  0  1 
0. 1764E+01 
0. 1 5  OOE  +  O 1 
0. 1220E+01 
0.9271 E+00 
0  «b  2  38E  +  00 
0.3136E+00 
0.1 036E-02 
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